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Executive  Summary 


This  report  presents  results  from  the  basic  contract  and  the  three  options.  The 
report  is  divided  into  2  parts.  This  approach  was  taken  to  organize  the  studies  and 
investigations  conducted  and  the  results  obtained  into  logical  reporting  categories.  In  the 
first  part,  we  first  review  the  approach  to  scene  partitioning  and  statistical  characterization 
originally  developed  at  Syracuse  University  by  Drs.  M.  A.  Slamani  and  D.  D.  Weiner. 
The  approach  consisted  previously  of  two  main  stages  referred  to  as  Mapping  and 
Indexing.  In  the  first  stage,  clutter  patches  are  separated  from  background  noise,  and,  in 
the  second  stage,  clutter  patches  are  separated  from  each  other  and  the  probability  density 
function  (PDF)  is  approximated  for  each  patch. 

The  approach  to  scene  partitioning  and  statistical  characterization  is  remodeled  so 
that  it  can  be  applied  in  general  to  separate  any  set  of  non-homogeneous  patches  without 
having  to  require  that  a  background  noise  be  present.  Also,  recognizing  that  two 
procedures  are  used  to  separate  regions,  the  stages  of  the  approach  have  been  redesigned 
to  reflect  this  fact.  The  two  main  stages  of  the  approach  are  known  as  (1)  the  Mapping 
procedure,  which  uses  image  processing  means  to  separate  regions,  and  (2)  the  Statistical 
procedure,  which  uses  statistical  means  to  separate  between  them.  The  approach  is 
assigned  the  acronym  A’ SCAPE  which  stands  for  the  Automated  Statistical 
Characterization  And  Partitioning  of  Environments.  Though  using  two  procedures, 
A’ SCAPE  needs  four  stages  to  achieve  its  aim  at  separating  the  different  contiguous  non- 
homogeneous  regions  that  may  exist  in  a  scene.  These  stages  consist  of  a  preprocessing 
stage  where  classical  time-frequency  processing  of  the  data  is  performed,  a  mapping 
procedure  stage,  a  statistical  procedure  stage,  and  an  indexing  stage  which  assigns  a  set  of 
descriptors  for  every  pixel  in  the  scene  under  investigation.  When  A’ SCAPE  is  followed 
by  a  detection  stage,  all  information  is  available  with  respect  to  which  probability  density 
function  (PDF)  should  be  used  for  each  cell  in  the  scene. 
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In  order  to  enable  A’ SCAPE  to  process  real  data,  the  image  processing  and 
statistical  algorithms  had  to  be  rewritten  and  new  rules  (in  addition  to  the  old  ones)  were 
selected  to  enable  the  convergence  of  the  expert  system  which  controls  the  flow  between 
the  different  stages  in  A’ SCAPE. 

An  example  of  a  real  IR  scene  is  included  in  the  report  which  consists  of  a  lake 
region  and  a  land  region.  First,  the  samples  are  selected  in  the  preprocessing  stage  to 
result  in  a  scene  with  uncorrelated  data.  Then,  the  mapping  procedure  separates  the  lake 
and  land  regions.  Road  ways  and  very  small  patches  in  the  land  region  are  successfully 
detected  as  outliers.  Also,  non-homogeneous  regions  (due  mainly  to  trees)  in  the  land 
region  near  the  boundary  between  the  land  and  lake  are  separated  by  the  statistical 
procedure  and  grouped  as  homogeneous  subpatches.  Finally,  every  homogeneous 
declared  region  has  its  PDF  approximated. 

A  demo  package  built  in  Matlab  was  also  developed  which  describes  in  detail  the 
different  stages  of  A’ SCAPE.  The  package  has  a  friendly  mouse  driven  graphical  user 
interface  (GUI)  and  consists  of  two  main  sections.  The  first  section  presents  the  detailed 
steps  of  the  real  IR  data  example.  The  second  section  is  subdivided  into  two  subsections 
where  the  first  one  consists  of  a  set  of  examples  which  illustrate  the  need  for  A’ SCAPE, 
whereas  in  the  second  a  description  is  given  for  every  stage  of  A’ SCAPE.  The  views  can 
be  displayed  manually  or  automatically.  A  movie  is  included  that  shows  the  steps  through 
which  the  scene  goes  during  the  mapping  procedure.  This  demo  package  is  among  the 
deliverables. 

The  objective  of  the  second  part  was  to  study  some  specific  parametric  methods 
for  the  detection  of  signals  in  the  presence  of  interferences  such  as  noise,  clutter  and 
jamming.  Specifically,  we  looked  at  a  technique  developed  by  LeCadre  and  studied  its 
performance  by  comparing  it  to  previously  developed  methods  and  against  several 
parameters. 

The  seope  of  this  effort  was  limited  to  the  feasibility  study  of  this  method  as  far  as 
its  implementation  is  concerned.  We  investigated  the  performance  of  this  technique  and 
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studied  all  other  effects  which  may  contribute  to  its  degradation  and  therefore  reduce  the 
probability  of  detection.  LeCadre’s  technique  is  based  upon  assuming  an  AR  or  (ARMA) 
model  of  the  interference.  The  first  step  in  the  technique  is  then  to  estimate  the  order  of 
the  model  used.  The  second  step  involves  estimating  the  AR  parameters  from  which  a 
covariance  matrix  is  derived.  Once  this  matrix  is  obtained,  it  will  be  used  in  conjunction 
with  the  detection  problem  to  study  two  direction  of  arrival  (DOA)  techniques,  namely 
MUSIC  and  ESPRIT,  using  the  transformed  or  sometimes  referred  to  as  the  “whitening” 
functional  obtained  from  the  inverse  of  the  AR  parameters  covariance  matrix. 

Note  that  in  this  effort,  we  extended  the  ideas  presented  by  LeCadre  to  include 
scenarios  where  signals,  clutter  and  additive  noise  are  included.  These  additions  make 
this  work  very  interesting  since  several  improvements  have  been  made  to  the  original 
work.  Moreover,  the  computer  simulations  of  section  3  show  that  the  results  do  match 
the  theoretical  development  of  section  2. 

We  feel  that  substantial  additions  to  the  beam-forming  problem  have  been  added 
in  this  effort  and  the  results  of  the  computer  simulations  describe  the  efficacy  of  the 
proposed  algorithms. 


Parti 

Automated  Statistical  Characterization  and  Partitioning  of  Environments 

with  Application  to  IR  Data 


1 


1 


INTRODUCTION 


1.1  -  Problem  Statement 

In  signal  processing  applications  it  is  common  to  assume  a  Gaussian  problem  in 
the  design  of  optimal  signal  processors.  However,  studies  have  shown  that  the  Gaussian 
reeeiver  performs  very  poorly  in  strong  interference  whenever  the  interference  and  signal 
spectra  carmot  be  separated  by  filtering.  For  example,  consider  the  spectra  shown  in 
Figure  1.1  consisting  of  24  Doppler  bins  with  uniformly  spaced  targets,  indicated  by  the 


Figure  1.1  -  Targets  in  Non-Gaussian  Interference 

small  arrows,  embedded  in  background  noise  and  a  bell  shaped  Gaussian  interference. 
The  optimal  Gaussian  based  detector,  referred  to  as  the  joint-domain  localized 
generalized  likelihood  ratio  receiver,  is  applied  to  each  Doppler  bin.  The  performance  of 


2 


the  receiver',  shown  in  Figure  1.2,  reveals  that  the  probability  of  detection,  PD,  of  the 
receiver  is  close  to  unity  everywhere  except  for  Doppler  bins  11,  12,  and  13  in  which  the 
strong  Gaussian  disturbance  exists  and  PD  falls  rapidly  to  the  probability  of  false  alarm, 
PFA. 


Because  of  the  fact  that  the  interference  is  Gaussian  and  that  the  detector  is  based 
on  the  Gaussian  distribution,  we  recognize  that  it  is  not  possible  to  obtain  any  detection 
improvement.  A  question  that  arises  is,  “Could  improved  detection  have  been  obtained  in 
bins  11  to  13  had  the  disturbance  been  non-Gaussian  ?”.  The  answer  is,  “It  is  possible  to 
achieve  significant  improvement  in  detection  performance  when  the  disturbance  is  non- 
Gaussian”.  To  illustrate  this  answer.  Table  1.1  presents  the  results  of  an  example  for  a 


'  H.  Wang,  L.  Cai,  "On  Adaptive  Multiband  Signal  Detection  with  GLR  algorithm."  IEEE  Trans,  on  AES, 
AES-27,  no.  2,  March  1991. 


3 


weak  target  in  strong  K-distributed  clutter  where  the  target  and  clutter  spectra  are 
coincident^  (i.e.,  not  separable).  Note  that  the  K-distributed  based  detector  provides 
significantly  improved  performance  even  when  the  PD  of  the  Gaussian  detector 
approaches  the  PFA.  For  example,  with  SCR=0  dB  and  PFA=10  ^  PD  equals  0.3  for  the 
k-distributed  detector  and  10"^  for  the  Gaussian  detector.  This  means  that  the  K- 
distributed  based  detector  can  detect  3  out  of  10  targets,  whereas  the  Gaussian  detector 
can  detect  only  1  in  1 00000  targets  on  the  average.  The  important  point  is  that  detections 
can  be  made  in  a  non-Gaussian  environment  using  a  non-Gaussian  detector  even  when 
such  performance  is  not  possible  with  the  Gaussian  detector.  Unfortunately,  there  are 
problems  associated  with  non-Gaussian  distributions:  (1)  there  are  a  multitude  of  Non- 
Gaussian  distributions,  (2)  some  non-Gaussian  distributions  have  shape  parameters  which 
result  in  different  shaped  probability  density  functions  (PDFs),  and  (3)  for  a  given  set  of 
data,  it  is  difficult  to  determine  which  PDF  can  approximate  the  data. 


SCR  (dB) 

PFA 

PD 

(K-Distributed 

Detector) 

PD 

(Gaussian 

Detector) 

0 

0.50 

0.06 

-10 

0.40 

0.02 

-20 

0.22 

0.01 

0 

10^ 

0.30 

lO'^* 

-10 

10-' 

0.25 

10-' 

-20 

10-' 

0.10 

10-' 

Table  1.1  -  Comparison  of  non-Gaussian  and  Gaussian  based  detectors 

The  use  of  PDF  based  detectors  in  the  implementation  of  likelihood  ratio  tests, 
LRTs,  and  locally  optimum  detectors,  LODs,  for  the  general  detection  problem  allows  us 


^  P.  Chakravathi,  M.  A.  Slamani,  D.  D.  Weiner,  "Performance  of  the  Locally  Optimum  Detector  in  a 
Correlated  K-Distributed  Disturbance,"  Proceedings  of  the  1993  National  Radar  Conference  on 
Revolutionary  Developments  in  Radar.  Wakefield,  MA,  April  20-23,  1993. 
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to  derive  algorithms  for  performing  both  strong  and  weak  signal  detection,  respectively, 
in  a  non-Gaussian  environment.  This  is  in  contrast  to  classical  detection  which  assumes  a 
priori  knowledge  of  the  joint  PDF  underlying  the  received  data  and  makes  use  of  a  single 
detector  usually  based  on  the  Gaussian  distribution. 

Note  that  in  practice,  one  cannot  assume  a  priori  knowledge  of  the  joint  PDF 
underlying  the  received  data  since  the  received  data  can  come  from  a  clear  region,  where 
background  noise  alone  is  present,  or  from  a  clutter  region,  where  returns  are  due  to 
reflections  from  such  objects  as  ground,  sea,  buildings,  birds,. ..etc.  When  a  desired  target 
return  is  from  a  clear  region  and  the  background  noise  is  sufficiently  small,  the  signal-to- 
noise  ratio  will  be  large  and  the  strong  signal  detector  (i.e,  the  LRT)  should  be  used. 
However,  if  a  desired  target  return  is  from  a  clutter  region,  two  situations  can  exist.  When 
the  desired  target  can  be  separated  from  the  clutter  by  means  of  space-time  processing 
and  the  background  noise  is  sufficiently  small,  the  signal-to-noise  ratio  will  be  large  and  a 
strong  signal  detector  should  again  be  used.  When  the  desired  target  cannot  be  separated 
from  the  clutter  by  means  of  space-time  processing  and  the  clutter  return  is  much  larger 
then  the  desired  target  return,  then  a  weak  signal  detector  (e.g.,  LOD)  based  on  the  PDF 
of  the  clutter  should  be  used.  Use  of  the  LOD  in  a  strong  signal  situation  can  result  in  a 
severe  loss  in  performance.  Hence,  it  is  necessary  for  the  receiver  to  determine  whether  a 
strong  or  weak  signal  situation  exists. 

Figure  1.3  summarizes  the  different  cases  which  may  arise  depending  on  whether 
the  target  is  embedded  in  background  noise  or  in  background  noise  plus  clutter.  This  can 
result  either  in  a  strong  signal  case,  intermediate  signal  case  or  a  weak  signal  case.  Note 
that  the  Gaussian  assumption  is  used  for  the  cases  where  the  likelihood-ratio-test  (LRT) 
and  generalized  likelihood-ratio-test  (GLRT)  detectors  are  utilized.  For  the  weak  signal 
case,  the  PDF  of  the  region  has  to  be  approximated  in  order  to  enable  the  use  of  the 
appropriate  LRT  or  the  sub-optimal  locally-optimum-detector  (LOD). 

All  of  this  suggests  the  necessity  for  a  procedure  to  1)  continuously  monitor  the 
environment,  2)  subdivide  the  surveillance  volume  into  homogeneous  patches,  and  3) 
select  the  appropriate  detector  for  processing  the  data.  In  addition,  depending  on 
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statistical  changes  in  the  environment  over  time  and  space,  the  process  enables  the 
receiver  to  adapt  so  as  to  obtain  close  to  optimal  performance.  This  is  achieved  by  the 
Automated  Statistical  Characterization  and  Partitioning  of  Environments  (A' SCAPE) 
procedure,  previously  used  successfully  on  simulated  radar  data^''. 
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Figure  1.3  -  Different  Target  Cases 


’  M.  A.  Slamani  and  D.  D.  Weiner,  "Rationale  Behind  the  Use  of  Image  Processing  to  Partition  a  Radar 
Surveillance  Volume  into  Background  Noise  and  Patches,"  Proceedings  of  the  1993  36th  Midwest 
Symposium  on  Circuits  and  Systems. 

“  M.  A.  Slamani  and  D.  D.  Weiner,  "Use  of  Image  Processing  to  Partition  a  Radar  Surveillance  Volume 
into  Background  Noise  and  Patches,"  Proceedings  of  the  1993  Conference  on  Information  Sciences  and 
Systems. 
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A' SCAPE  uses  two  separate  procedures  to  determine  all  homogeneous  regions 
and  sub-regions  in  the  scene.  The  first  procedure,  referred  to  as  the  mapping  procedure,  is 
used  to  separate  contiguous  homogeneous  regions  by  segregating  between  their  power 
levels.  The  second  procedure,  referred  to  as  the  statistical  procedure,  separates  contiguous 
homogeneous  regions  by  segregating  between  their  data  distributions.  The  statistical 
procedure  uses  the  Ozturk  algorithm,  a  newly  developed  algorithm  for  analyzing  random 
data^  Furthermore,  the  statistical  procedure  identifies  suitable  approximations  to  the 
probability  density  function  for  each  region.  Convergence  of  the  mapping  and  statistical 
procedures  are  controlled  through  expert  system  rules  as  developed  in  Chapter  4. 

1.2  -  A’SCAPE  Procedure 

A’ SCAPE  is  based  on  the  approach  to  scene  partitioning  and  statistical 
characterization  originally  developed  at  Syracuse  University  by  Drs.  M.  Adel  Slamani 
and  D.  D.  Weiner.  The  approach  consisted  previously  of  two  main  stages  referred  to  as 
Mapping  and  Indexing.  In  the  first  stage,  clutter  patches  are  separated  from  background 
noise,  and,  in  the  second  stage,  contiguous  non-homogeneous  clutter  patches  are 
separated  from  each  other  and  the  probability  density  function  (PDF)  is  approximated  for 
each  patch. 

The  approach  is  remodeled  so  that  it  can  be  applied  in  general  to  separate  any  set 
of  contiguous  non-homogeneous  patches  without  having  to  require  that  a  background 
noise  be  present.  Also,  recognizing  that  two  procedures  are  used  to  separate  regions,  the 
stages  of  the  approach  have  been  redesigned  to  reflect  this  fact.  The  two  main  stages  of 
the  approach  are  known  as  (1)  the  Mapping  procedure,  which  uses  image  processing 
means  to  separate  between  regions,  and  (2)  the  Statistical  procedure,  which  uses 
statistical  means  to  separate  between  regions.  The  approach  is  assigned  the  acronym 
A’ SCAPE  which  stands  for  the  Automated  Statistical  Characterization  And  Partitioning 
of  Environments.  Though  using  two  procedures.  A’ SCAPE  needs  four  blocks  to  achieve 
its  goal  of  separating  the  different  contiguous  non-homogeneous  regions  that  may  exist  in 

^  A.  Ozturk  and  E.  J.  Dudewicz,  "A  New  Statistical  Goodness  of  Fit  Test  Based  on  Graphical 
Representation,"  The  Biomedical  Journal,  1991. 
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a  scene.  As  shown  in  Figure  1.4,  the  first  block  is  a  preprocessing  block  that  performs 
classical  space-time  processing  on  the  collected  data  (e.g.,  inphase  and  quadrature 
components  extraction,  envelope  detection).  The  second  block  separates  contiguous 
homogeneous  patches  and  subpatches  by  segregating  between  their  average  power  levels 
based  on  the  mapping  procedure  described  in  chapter  2.  The  next  block,  presented  in 
chapter  3,  goes  one  step  further  and  separates  homogeneous  subpatches  by  segregating 
between  their  probabilistic  data  distributions.  Furthermore,  this  block  identifies  suitable 
approximations  to  the  probability  density  function  (PDF)  of  each  homogeneous  patch  and 
determines  the  location  of  outliers  in  the  scene.  The  final  block,  labeled  indexing,  indexes 
the  scene  under  investigation  by  assigning  a  set  of  descriptors  to  every  pixel  in  the  scene. 
For  each  pixel,  the  indexing  is  used  to  indicate  to  which  homogeneous  patch  the  pixel 
belongs,  whether  it  is  an  edge  pixel,  whether  it  is  an  outlier,  which  pixels  can  be  chosen 
as  reference  pixels  if  the  pixel  is  to  be  tested,  and  which  PDF  best  approximates  the  data 
in  the  pixel.  Note  that  the  reference  pixels  are  pixels  which  belong  to  the  same 
homogeneous  patch  and  which  are  closest  to  the  pixel  to  be  tested. 


Preprocessing 


Patch  &  Subpatch 
Investigation  Using  the 
Mapping  Procedure 

t 

Outlier  Search  & 
Subpatch  Investigation 
Using  the 

Statistical  Procedure 


Indexing 


Figure  1.4  -  A’SCAPE  Block  Diagram 
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The  forward  and  backward  interactions  between  the  different  blocks  are 
controlled  through  rules  of  an  expert  system  shell  referred  to  as  Integrated  Processing  and 
Understanding  of  Signals  (IPUS)  developed  jointly  by  the  University  of  Massachusetts 
and  Boston  University^  IPUS  is  discussed  in  chapter  4. 

Note  that  when  A'SCAPE  is  followed  by  a  detection  stage  (e.g.,  target  detection  in 
Radar),  all  information  needed  is  available  for  every  pixel  in  the  scene.  Furthermore, 
given  the  PDF  that  can  approximate  the  patch  in  which  the  test  pixel  is  located,  the 
appropriate  detector  is  readily  selected.  This  is  in  contrast  to  the  classical  detection 
approach  where  a  single  detector  (usually  the  Gaussian  detector)  is  used  in  processing  the 
entire  scene. 

The  mapping  and  statistical  procedures  are  presented  in  chapters  2  and  3, 
respectively.  The  expert  system  shell  is  discussed  in  chapter  4.  An  example  illustrating 
the  different  stages  of  A'SCAPE  when  applied  to  real  data  of  an  IR  image  is  given  in 
chapter  5. 


^  H.  Nawab,  V.  Lesser,  et  al.,  Integrated  Processing  and  Understanding  of  Signals  in  Symbolic  and 
Knowledge-Based  Signal  Processing.  Edited  by  A.V.  Oppenheim  and  H.  Nawab,  Prentice  Hal, 
Ennglewood  RPsiffs,  N.J.,  1991. 
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2 


Mapping  Procedure 


2.1  -  Introduction 

In  this  chapter,  a  mapping  procedure  is  presented  to  partition  a  scene  into 
homogeneous  regions  by  segregating  between  the  power  levels  of  the  different  patches. 

Given  a  scene  which  consists  of  a  set  of  non-homogeneous  patches,  the  mapping 
procedure  starts  first  by  separating  the  regions  into  a  patch  with  the  lowest  power  level, 
referred  to  as  the  lowest  patch  (LP),  and  remaining  patches  (RPs).  The  procedure  is  then 
repeated  to  isolate  the  next  patch  with  the  lowest  power  level  (LP)  and  so  on  until  it  is  not 
possible  to  find  anymore  LPs.  General  observations  are  first  made  on  LP  and  RPs. 

2.1.1  -  Observations  on  LP  and  RPs 

Assume  that  a  scene  consists  of  JxK  pixels  and  that  JxK  data  magnitudes  PG,k) 
are  available  to  the  mapping  procedure.  In  this  case,  JxK  pixels  need  to  be  mapped  into 
LP  and  RPs.  Let's,  examine  the  nature  of  the  LP  and  RPs  in  order  to  understand  the 
theory  behind  the  procedure  developed  for  mapping. 

The  following  observations  are  based  on  computer  generated  examples  of  LP  and 
RPs  data  where  the  RPs-to-LP  (RLR)  ratio  is  assumed  to  be  greater  than  0  dB. 

2.1.1. a  -  Observations  on  LP 

-  On  average,  the  LP  data  values  are  smaller  than  the  RPs  data  values. 

-  Large  data  values  exist  in  a  LP  that  may  be  higher  than  some  data  values  of  the 

RPs. 

-  Large  data  values  in  the  LP  tend  to  be  isolated  points. 

-  The  number  of  LP  data  significantly  larger  than  the  average  is  relatively  small. 
Figure  2.1  shows  a  typical  LP  data  histogram. 

-  The  relatively  small  number  of  LP  pixels  with  large  data  values  are  distributed 
evenly  throughout  the  LP. 
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Frequency 


Figure  2.1  -  Typical  LP  data  histogram 

2.1. l.b  -  Observations  on  RPs 

-  On  average,  RPs  data  values  are  higher  than  LP  data  values. 

-  The  large  RPs  data  values  are  larger  than  the  largest  LP  data  values  assuming 
positive  RPs-to-LP-Ratio  (RLR). 

-  Small  RPs  data  values  exist  and  may  be  smaller  than  the  large  LP  data  values. 

-  Large  data  values  in  the  RPs  tend  to  be  clustered. 

2.2  -  Mapping  Procedure 

Using  the  fact  that  the  RPs  patches,  on  average,  have  stronger  magnitudes,  the 
mapping  procedure  begins  by  setting  a  threshold  that  results  in  a  specified  fraction  of  LP 
pixels.  Image  processing  tools  based  on  data  thresholding  and  region  partitioning  are  then 
used  to  establish  the  LP  and  RPs  patches.  If  the  final  image  contains  a  significantly 
different  fraction  of  LP  than  originally  established  by  the  initial  threshold,  the  process  is 
repeated  with  a  new  threshold.  The  mapping  procedure  iterates  until  it  is  satisfied  that  the 
final  scene  is  consistent  with  the  latest  specified  threshold.  Then,  edges  of  all  patches  are 
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detected  using  an  image  processing  technique  based  on  the  unsharp  masking’^  Finally,  all 
regions  in  the  scene  are  assigned  a  unique  identifying  number. 

The  mapping  procedure  is  therefore  composed  of  three  steps.  The  first  step  results 
on  the  identification  of  LP  and  RPs,  whereas  the  second  and  third  steps  consist  of  the 
enhancement  and  detection  of  patch  edges  and  labeling  of  the  different  regions, 
respectively.  These  three  steps  are  explained  next. 

2.2.1  -  Identification  of  LP  and  RPs 

Identification  of  LP  and  RPs  is  performed  by  the  following  steps:  thresholding, 
quantization,  corrections  and  assessment. 

2.2.1. a  -  Thresholding  and  Quantization 

Identification  of  LP  and  RPs  starts  by  setting  a  threshold  q  that  results  in  a 
specified  fraction  of  LP  declared  pixels.  Then,  a  quantized  volume  is  formed  as  follows: 
all  pixels  with  magnitude  less  than  q  are  given  a  value  of  0  (zero)  and  all  pixels  with 
magnitude  above  q  are  given  a  value  of  1  (one).  Let  Q(j,k)  represent  the  quantized  value 
of  the  jkth  pixel.  Then, 


f  1  if  Pij,k)  >  q 

n  /P  ^  rl,2,...,Jandk=l,2,...,K  (2.1) 

[0  if  P{j,k)<q 

where  P(j,k)  is  the  magnitude  of  the  jkth  pixel. 

2.2.1.b  -  Corrections 

Consider  a  set  of  3x3  pixels.  As  shown  in  Figure  2.2,  let  the  center  pixel  be 
referred  to  as  the  test  pixel  and  the  surrounding  pixels  be  referred  to  as  the  neighboring 

’  R.  Gonzalez,  P.  Wintz,  "Digital  Image  Processing."  2nd  edition,  Addison- Wisley  Publishing  Company, 
Nov.  1987. 

*  E.  Hall,  Computer  Image  Processing  and  Recognition.  Academic  Press,  1979. 
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pixels.  Assume  that  a  patch  cannot  be  formed  by  a  single  pixel.  In  this  case,  every  test 
pixel  in  a  patch  has  at  least  one  neighboring  pixel  that  belongs  to  the  same  patch. 

A  test  pixel  belonging  to  a  RPs  patch  that  has  at  least  one  neighboring  LP  pixel  is 
referred  to  as  a  RPs  edge  pixel  (RPsE).  On  the  other  hand,  a  test  pixel  that  belongs  to  a 
RPs  patch  for  which  none  of  the  neighboring  pixels  are  in  the  LP  is  referred  to  as  an  inner 
RPs  pixel.  Vice  versa,  a  test  pixel  belonging  to  a  LP  patch  that  has  at  least  one 
neighboring  RPs  pixel  is  referred  to  as  a  LP  edge  pixel  (LPE).  Also,  a  test  pixel  that 
belongs  to  a  LP  patch  for  which  none  of  the  neighboring  pixels  are  in  the  RPs  is  referred 
to  as  an  iimer  LP  pixel. 


Figure  2.2  -  3x3  pixels 

The  proposed  correction  technique  consists  of  transforming  the  quantized  volume 
into  a  "corrected"  volume.  The  transformation  consists  of  the  following  steps: 

- 1.  In  the  quantized  volume,  declare  as  LP  pixels  all  pixels  with  quantized  values 
Q(j,k)=0  and  as  RPs  pixels  all  pixels  with  quantized  values  Q(j,k)=l. 
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-  2.  Choose  the  necessary  number  of  RPs  neighboring  pixels,  NCQ,  for  a  test  pixel 
in  the  quantized  volume  to  be  declared  as  a  RPs  pixel  in  the  corrected  volume.  NCQ  can 
take  one  of  the  following  values:  5,6, 7,8. 

-  3.  For  every  test  pixel  in  the  quantized  volume  count  the  number  of  neighboring 
RPs  pixels.  If  the  number  is  greater  than  or  equal  to  NCQ  declare  the  test  pixel  as  a  RPs 
pixel  in  the  corrected  volume.  Otherwise,  declare  the  test  pixel  as  a  LP  pixel  in  the 
corrected  volume. 

When  all  the  pixels  of  the  quantized  volume  have  been  tested,  a  "corrected" 
volume  consisting  of  declared  LP  or  RPs  pixels  is  obtained. 

Because  NCQ  is  chosen  to  be  relatively  large  (i.e.  NCQ=5,6,7  or  8),  LP  pixels 
that  were  incorrectly  identified  in  the  quantized  volume  as  RPs  pixels  due  to  their  large 
power  tend  to  be  reclassified  as  LP  pixels.  Also,  inner  RPs  pixels  in  the  quantized  volume 
are  recognized  as  RPs  pixels  in  the  "corrected"  volume.  Meanwhile,  most  of  the  RPs  edge 
pixels  in  the  quantized  volume  are  recognized  as  LP  pixels  in  the  "corrected"  volume. 
This  results  in  an  over-correction  where  most  of  the  RPs  edge  pixels  are  identified  as  LP 
pixels.  As  an  example,  when  NCQ=8,  only  inner  RPs  pixels  in  the  quantized  volume  are 
recognized  as  RPs  pixels  in  the  "corrected"  volume  and  all  RPs  edge  pixels  in  the 
quantized  volume  are  recognized  as  LP  pixels  in  the  "corrected"  volume.  In  order  to 
recover  the  edge  pixels,  a  second  correction  stage  is  needed  where  the  "first  corrected" 
volume  will  be  transformed  into  a  "second  corrected"  volume.  Let  the  "first  corrected" 
volume  be  referred  to  as  the  "corrected-quantized"  volume  (CQV)  and  the  "second 
corrected"  volume  be  referred  to  as  the  "corrected-corrected"  volume  (CCV).  The 
following  steps  are  used  to  transform  the  CQV  into  the  CCV : 

- 1.  Choose  the  necessary  number  of  RPs  neighboring  pixels,  NCC,  for  a  test  pixel 
in  the  CQV  to  be  declared  as  a  RPs  pixel  in  the  CCV.  NCC  can  take  one  of  the  following 
values:  1,2,3  or  4. 

-  2.  For  every  test  pixel  in  the  CQV  count  the  number  of  neighboring  RPs  pixels. 
If  the  number  is  greater  than  or  equal  to  NCC,  declare  the  test  pixel  as  a  RPs  pixel  in  the 
CCV.  Otherwise,  declare  the  test  pixel  as  a  LP  in  the  CCV. 
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2.2.1.C  -  Assessment 


Let  LPQP,  LPCQP  and  LPCCP  denote  the  percentage  of  LP  pixels  in  the 
quantized,  "corrected-quantized"  and  "corrected-corrected"  volumes,  respectively.  LPQP 
is  pre-specified  so  as  to  determine  the  threshold  for  the  quantized  volume,  whereas 
LPCQP  and  LPCCP  are  computed  after  the  CQV  and  the  CCV  are  obtained. 

The  assessment  process  consists  of  comparing  LPCQP  and  LPCCP  to  LPQP  in 
order  to  determine  whether  or  not  the  percentages  of  the  LP  pixels  after  correction  are 
consistent  with  the  percentage  of  LP  pixels  in  the  quantized  volume.  When  there  is  no 
consistency,  further  quantization,  correction  and  assessment  are  performed  until 
consistency  is  obtained. 

The  thresholding/quantization,  first  correction,  second  correction,  and  assessment 
stages  are  used  to  find  the  best  threshold  to  separate  between  LP  and  RPs  patches.  Once 
LPQP  has  been  set  a  threshold  is  computed.  Then,  corrections  are  made  to  try  and  build 
the  LP  region  and  the  RPs  patches  .  The  correction  stages  re-label  some  of  the  above¬ 
threshold  pixels  as  LP  pixels  if  they  are  likely  to  belong  to  the  LP,  and  some  of  the 
below-threshold  pixels  as  RPs  pixels  if  they  are  likely  to  belong  to  a  RPs  patch,  based  on 
the  ehoices  for  NCQ  and  NCC.  Depending  on  the  quality  of  the  threshold  choice  many  or 
few  pixels  are  re-labeled.  At  the  end  of  the  procedure,  LPCCP  is  computed  and  compared 
to  LPQP,  if  the  values  are  within  a  certain  range,  few  pixels  would  have  been  re-labeled, 
the  threshold  is  accepted  and  the  assessment  passes.  Otherwise,  many  pixels  would  have 
been  re-labeled  and  the  threshold  is  rejected.  The  iterative  process  continues  then  by 
setting  another  threshold  through  the  choice  of  a  new  value  for  LPQP. 

Rules  for  choosing  NCQ,  NCC  and  LPQP  and  for  determining  when  consistency 
of  the  percentages  is  obtained  are  explained  in  Chapter  4. 

2.2.2  -  Enhancement  and  Detection  of  Patch  Edges 

Once  the  LP  and  RPs  have  been  detected.  The  mapping  procedure  enhances  the 
patches  and  detects  their  edges.  These  are  done  by  the  smoothing,  edge  detection  and 
edge  enhancement  stages  presented  next. 
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2.2.2.a  -  Smoothing 


In  many  cases  of  simulated  data,  examples  have  shown  that  when  the  percentages 
are  consistent,  declared  patches  may  contain  isolated  LP  declared  pixels.  Because  small 
magnitudes  can  arise  in  a  RPs  patch  as  explained  in  the  Section  2.1.1,  it  is  most  likely 
that  the  LP  isolated  pixels  in  the  RPs  patches  are  RPs  pixels.  The  smoothing  process  is 
used  to  detect  these  isolated  pixels  and  label  them  adequately  by  transforming  the  CCV 
into  a  smoothed  volume  (SV).  The  smoothing  technique  consists  of  the  following  steps: 

-  1.  Choose  the  necessary  number  of  RPs  neighboring  pixels  NS  for  a  LP 
identified  test  pixel  in  the  CCV  to  be  declared  as  a  RPs  pixel  in  the  SV  where  NS  can 
take  one  of  the  following  values:  5,6,7,  or  8. 

-  2.  For  every  LP  identified  pixel  in  the  CCV  count  the  number  of  neighboring 
RPs  pixels.  If  the  number  is  greater  than  or  equal  to  NS,  declare  the  test  pixel  as  a  RPs 
pixel  in  the  SV.  Otherwise  declare  the  test  pixel  as  a  LP  pixel  in  the  SV. 

2.2.2.b  -  Detection  of  Patch  Edges  and  Edge  Enhancement 

(i)  -  Detection  of  Patch  Edges 

After  smoothing,  each  pixel  in  the  SV  has  been  declared  as  either  a  LP  or  a  RPs 
pixel.  The  next  step  is  to  determine  which  of  the  RPs  pixels  are  located  on  the  edges  of 
the  RPs  patches  and  which  of  the  LP  pixels  are  located  on  the  edges  of  the  LP  patch.  This 
is  important  for  subsequent  signal  processing  if  reference  pixels  for  estimating  parameters 
of  a  test  pixel  are  to  be  chosen  properly. 

Identification  of  RPs  edge  pixels  (RPsE)  and  LP  edge  pixels  (LPE)  is  done  by  the 
use  of  an  image  processing  technique  referred  to  in  the  image  processing  literature  as  the 
unsharp  masking.  It  consists  of  the  following  steps: 

To  detect  RPsE  pixels 

1  -  A  weighting  filter  consisting  of  a  3x3  array  of  pixels  is  constructed,  as  shown 
in  Figure  2.3,  where  the  center  pixel  has  a  weight  given  by  w(0,0)=8  and  the  neighboring 
pixels  have  weights  given  by  w(-l,-l)=w(0,-l)=w(l,-l)=w(-l,0)=w(l,0)= 
w(-l,l)=w(0,l)=w(l,l)=-L  The  center  pixel  is  positioned  on  the  test  pixel.  Notice  that 
the  weights  of  the  filter  pixels  sum  to  zero.  In  particular. 
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2  =  0  (2.2) 

m=~]  n=-] 

2  -  Assume  that  the  weighting  filter  is  centered  at  the  jkth  pixel  in  SV.  The  pixels 
corresponding  to  the  3x3  array  of  the  weighting  filter  have  quantized  values  as  illustrated 
in  Figure  2.4.  By  definition. 


SQ(j,k) 


fl  if  the  jkth  cell  inSV  is  declared  as  RPs 
[  0  if  the  jkth  cell  in  SV  is  declared  as  LP 


(2.3) 


where  j=l,2,...,J  and  k=l,2,...,K. 

To  avoid  filter  pixels  falling  outside  SV,  the  coordinates  of  the  jkth  pixel  at  which 
the  filter  is  centered  are  constrained  to  j=2,3,...,J-l,  and  k=2,3,...,K-l. 

3  -  Form  the  sum, 


1  1 

w(m,n)  SQ(j  +  m,k  +  n)  (2.4) 

Wl=-1  W=-l 

-  If  S  is  equal  to  zero,  all  pixels  have  the  same  assigned  value.  This  can 
arise  only  when  the  test  pixel  is  not  an  edge  pixel. 

-  If  S  is  positive,  the  test  pixel  is  an  edge  pixel  and  is  labeled  as  such. 

-  If  S  is  negative,  the  test  pixel  cannot  be  an  edge  pixel.  On  the  other  hand, 
one  or  more  of  the  neighboring  pixels  are  guaranteed  to  be  an  edge  pixel. 

The  three  situations  are  illustrated  in  Figures  2.5-a,  b,  c  and  d.  In  Figures  2.5-a  and  2.5-b, 
S=0  because  all  9  pixels  are  in  LP  and  RPs,  respectively.  Observe  that  the  test  pixel  is  not 
an  edge  pixel.  In  Figure  2.5-c,  S=4>0.  Note  that  the  test  pixel  is  an  edge  pixel.  Finally,  in 
Figure  2.5-d,  S=-2<0  and  the  test  pixel  is  not  an  edge  pixel. 


To  detect  LPE  pixels. 


First,  Equation  2.3  is  rewritten  as  follows. 


SQU,k) 


if  the  jkth  cell  in  SV  is  declared  as  RPs 
|l  z/  the  jkthcell inSVisdeclared as  LP 


where  j=l,2,...,J  and  k=l,2,...,K. 


(2.5) 
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Figure  2.5  -  Example  of  Unsharp  Masking  Technique 
(a)  S=0,  (b)  S=0,  (c)  S>0,  (d)  S<0 

Then,  step  3  above  is  repeated. 

At  the  end  of  the  edge  detection  procedure,  each  pixel  in  the  original  volume  has 
been  labeled  as  either  RPs,  LP,  RPsE,  or  LPE. 


(ii)  Enhancement  of  Patch  Edges 

The  edges  of  the  smoothed  volume  (SV)  tend  not  to  follow  the  irregular  edges 
that  may  actually  exist.  Consequently,  the  edges  are  further  enhanced  by  examining  the 
magnitudes  of  pixels  just  outside  the  LPE  pixels  and  the  RPsE  pixels.  If  the  magnitudes 
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of  these  pixels  exceed  the  last  threshold  q  for  which  the  assessment  passes,  they  are 
declared  as  RPsE  pixels  otherwise  they  are  declared  as  LPE  pixels. 

At  the  end  of  the  edge  enhancement  procedure,  edges  are  detected  and  each  pixel 
in  the  original  volume  is  labeled  as  either  LP,  RPs,  LPE,  or  RPsE  pixel. 

2.2.3  -  Labeling  of  Patches 

When  the  first  round  of  the  mapping  procedure  is  completed,  recall  that  the 
mapped  volume  has  a  value  of  0  assigned  to  LP  pixels  and  a  value  of  1  assigned  to  the 
RPs  pixels.  Therefore,  nothing  more  needs  to  be  done  for  the  LP  region  as  all  of  its  pixels 
are  already  indexed  by  the  number  0.  On  the  other  hand,  all  pixels  in  each  of  the  RPs 
patches  are  assigned  a  value  of  1 .  Thus,  a  numbering  procedure  has  to  be  implemented  to 
enable  the  computer  to  distinguish  between  the  various  RPs  patches.  The  approach  taken 
in  this  work  is  to  assign  every  pixel  in  the  first  patch  investigated  the  number  1,  every 
pixel  in  the  second  patch  investigated  the  number  2,  and  so  on  until  all  patches  have  been 
indexed  with  consecutive  integers.  In  this  way,  all  pixels  in  each  RPs  patch  are  assigned  a 
unique  number. 

If  a  pixel  belongs  to  a  new  RPs  patch,  the  key  to  the  numbering  is  the  ability  to 
recognize  this  fact.  This  is  done  by  considering  a  mask  of  5  pixels  as  shown  in  Figure  2.6 
where  the  white  pixels  represent  neighboring  pixels  and  the  shaded  one  is  the  test  pixel  to 
be  numbered. 

If  M(j,k)  is  the  value  assigned  to  the  ijth  pixel  in  the  mapped  volume,  then 


M(j,k)  = 


0  if  the  jkth  pixel  is  declared  as  LP 
1  if  the  jkth  pixel  is  declared  as  RPs 


(2.6) 


Assuming  that  the  test  pixel  to  be  numbered  is  the  jkth  pixel  in  the  original  surveillance 
volume,  let  the  assigned  number  be  denoted  by  N(j,k).  Also,  let  G  denote  the  unique 
number  assigned  to  the  RPs  patch  previously  investigated  and  H  the  minimum  positive 
number  assigned  to  the  neighboring  pixels.  Then,  by  definition,  we  have  that. 
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Figure  2.6  -  Mask  used  in  numbering 


NU,k) 


0  IfMU,k)  =  Q 

■G  +  1  If  all  neighboring  cells  are  numbered  0  (2.6) 

H  If  at  least  one  of  the  neighboring  cells  is  numbered  nonzero 


The  number  G  is  incremented  by  unity  whenever  a  new  RPs  patch  is  detected. 

Because  a  RPs  patch  boundary  may  be  sharply  shaped,  as  shown  in  the  example 
of  Figure  2.7,  the  numbering  procedure  may  end  up  by  assigning  two  different  numbers 
for  different  pixels  of  the  same  RPs  patch.  This  anomaly  is  avoided  by  further  testing  the 
neighboring  pixels  of  the  pixel  to  be  numbered  as  follows: 

1  -  For  a  given  pixel  to  be  numbered,  look  up  the  numbers  assigned  to  the  set  of 

neighboring  pixels  (j-l,k-l),  (j+l^k-l),  and  G'kk), 

2  -  Take  the  minimum  nonzero  number  of  those  in  step  1, 

3  -  Reassign  all  nonzero  numbered  neighboring  pixels  the  minimum  nonzero 
number  from  step  2, 
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4  -  Revisit  all  the  pixels  in  the  surveillance  volume  that  have  been  numbered  so 
far.  If  any  pixel  is  assigned  a  nonzero  number  identical  to  one  of  those  in  step  1,  reassign 
that  pixel  the  minimum  nonzero  number  of  step  2. 

For  example,  with  respect  to  Figure  2.7,  the  above  steps  have  the  effect  of 
assigning  a  value  of  1  to  all  pixels  of  the  RPs  patch  shown. 

Once  numbering  is  completed,  the  LP  pixels  are  assigned  a  value  of  0,  and  every 
RPs  patch  is  assigned  a  unique  positive  number. 


Figure  2.7  -  Example  of  a  sharply  shaped  boundary 


2.3  -  Conclusion 

The  mapping  procedure  consists  of  the  following  steps:  thresholding/quantization, 
correction,  assessment,  smoothing,  edge  enhancement,  edge  detection,  and  labeling  of 
patches.  As  explained  in  the  previous  sections  and  shown  in  Figure  2.8,  these  are 
subdivided  into  three  main  stages  referred  to  as  the  Identification  of  LP  and  RPs,  the 
detection  of  patch  edges,  and  labeling  of  patches,  respectively. 
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Identification  of 
LP  and  RPs 


Detection  of 
Patch  Edges 


Labeling  of 
Patches 


Figure  2.8  -  Block  Diagram  of  the  Mapping  Procedure 


As  stated  before,  once  numbering  is  completed,  the  LP  pixels  are  assigned  a  value 
of  0,  and  every  RPs  patch  is  assigned  a  unique  positive  number. 

Recall  that  mapping  consists  of  appropriately  selecting  a  threshold  to  distinguish 
between  LP  and  RPs  patches  using  only  the  assumption  that  the  LP  magnitudes,  on 
average,  are  smaller  than  the  RPs  magnitudes.  This  same  approach  may  be  used  once 
again  to  extract  a  LP  subpatch  from  a  set  of  contiguous  RPs  subpatches  of  higher 
magnitudes  in  a  given  RPs  patch.  In  this  case,  the  RPs  patch  containing  non- 
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homogeneous  RPs  subpatches  will  be  viewed  as  a  volume  containing  a  RPs  subpatch 
region  with  low  power  level  and  a  set  of  subpatches  with  higher  power  levels  occupying 
the  rest  of  the  RPs  patch  area. 

The  mapping  procedure  is  therefore  used  to  extract  the  RPs  subpatch  with  the 
lowest  power  (which  becomes  an  LP)  from  among  the  remaining  RPs  subpatches  in  a 
given  RPs  patch.  Because  the  numbering  stage  has  already  labeled  each  patch  with  a 
unique  number  it  is  straight  forward  for  the  program  to  select  a  patch  and  check  for  the 
presence  of  subpatches  in  it. 

For  each  RPs  patch,  the  mapping  procedure  is  performed  iteratively  until  it  is 
hypothesized  that  every  subpatch  in  a  given  RPs  patch  is  homogeneous  and  cannot  be 
partitioned  further.  After  all  RPs  subpatches  have  been  extracted,  the  surveillance  volume 
will  consist  of  a  set  of  homogeneous  patches  with  different  power  levels. 
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3 


Statistical  Procedure 


3.1  -  Introduction 

In  Chapter  2,  a  mapping  procedure  was  presented  to  partition  a  scene  into 
homogeneous  regions  by  segregating  between  the  power  levels  of  the  different  patches. 
When  no  more  patches  (sub-patches)  can  be  found,  the  mapping  procedure  ends  and  is 
followed  by  the  statistical  procedure.  The  statistical  procedure  is  applied  to  every 
homogeneous  previously  declared  patch  or  subpatch  to  (1)  further  separate  non- 
homogeneous  subpatches  having  very  similar  power  levels  but  different  statistical 
distributions,  (2)  locate  outliers  in  the  scene,  and  (3)  approximate  the  probability  density 
function  (PDF)  of  each  homogeneous  region. 

The  Ozturk  algorithm  is  used  by  the  statistical  procedure  to  approximate  the  PDF 
of  each  patch  and  is  presented  next  followed  by  the  definition  of  outliers.  The  strategy 
used  in  the  statistical  procedure  is  discussed  in  the  last  subsection. 

3.2  -  Introduction  to  the  Ozturk  Algorithm 

The  Ozturk  algorithm  is  based  on  sample  order  statistics  and  is  used  for  univariate 
distribution  approximation.  This  algorithm  has  two  modes  of  operation.  In  the  first  mode, 
the  algorithm  performs  a  goodness  of  fit  test.  The  test  determines,  to  a  desired  confidence 
level,  whether  the  random  data  is  statistically  consistent  with  a  specified  probability 
distribution.  In  the  second  mode  of  operation,  the  algorithm  approximates  the  PDF 
underlying  the  random  data.  In  particular,  by  analyzing  the  random  data  and  without  any 
a  priori  knowledge,  the  algorithm  identifies  from  a  stored  library  of  PDFs  that  density 
function  which  best  approximates  the  data.  Estimates  of  the  location,  scale,  and  shape 
parameters  of  the  PDF  are  provided  by  the  algorithm.  The  algorithm  has  been  found  to 
work  reasonably  well  for  observation  sizes  as  small  as  100.  Throughout  this  work,  100 
reference  pixels  are  selected  for  approximating  the  PDF  of  each  test  pixel.  Note  that  when 
a  region  contains  less  than  100  pixels  it  is  discarded  and  is  not  processed  by  the  statistical 
procedure  stage. 
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3.2.1  -  Goodness  of  Fit  Test 


The  goodness  of  fit  test  determines  whether  or  not  the  set  of  data  samples 
provided  to  the  algorithm  is  statistically  consistent  with  a  specified  distribution,  referred 
to  as  the  null  hypothesis.  Let  NR  denote  the  number  of  reference  pixels.  For  the  null 
hypothesis,  the  program  utilizes  a  Monte-Carlo  simulation  of  2,000  trials  to  generate  an 
averaged  set  of  NR  linked  vectors  in  the  UV  plane,  as  shown  in  Figure  3.1.  Using  the 
standardized  sample  order  statistics  of  the  data,  the  program  then  creates  a  second  system 
of  NR  linked  vectors  in  the  UV  plane.  The  terminal  points  of  the  linked  vectors,  as  well 
as  the  shapes  of  their  trajectories,  are  used  in  determining  whether  or  not  to  accept  the 
null  hypothesis.  The  null  hypothesis  is  the  distribution  against  which  the  sample  data  is  to 
be  tested. 


Figure  3.1  -  Goodness  of  fit  test 

The  algorithm  provides  quantitative  information  as  to  how  consistent  the  sample 
data  set  is  with  the  null  hypothesis  distribution  by  use  of  confidence  contours  where  each 
contour  is  derived  from  a  specified  probability  that  the  end  point  falls  within  the  contour 
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given  that  the  data  comes  from  the  null  distribution.  An  example  of  these  contours  is 
shown  in  Figure  3.1  for  probabilities  of  0.9,  0.95,  and  0.99.  If  the  end  point  of  the  sample 
data  linked  vector  locus  falls  within  a  contour,  then  the  sample  data  set  is  said  to  be 
statistically  consistent  with  the  null  hypothesis  at  a  confidence  level  based  on  the 
probability  specified  for  that  contour.  If  the  sample  data  set  is  truly  consistent  with  the 
null  hypothesis,  the  system  of  sample  linked  vectors  is  likely  to  closely  follow  that  for  the 
system  of  null  linked  vectors. 

3.2.2  -  Approximation  Chart  Mode 

The  approximation  chart  mode  is  simply  an  extension  of  the  goodness  of  fit  test 
mode.  Following  a  similar  approach  to  that  outlined  in  the  section  for  the  goodness  of  fit 
mode,  random  samples  are  generated  from  a  library  of  different  univariate  probability 
distributions.  In  the  goodness  of  fit  test  mode,  the  locus  end  point  was  obtained  for  the 
null  hypothesis  and  sample  size,  NR.  For  the  approximation  chart  mode  we  go  one  step 
further  by  obtaining  the  locus  end  point  for  each  distribution  from  the  library  of 
distributions  for  the  given  sample  size,  NR,  and  for  various  choices  of  the  shape 
parameter(s).  Thus,  depending  on  whether  it  has  a  shape  parameter  or  not,  each 
distribution  is  represented  by  a  trajectory  or  point  in  the  two  dimensional  UV  plane. 
Figure  3.2  shows  an  example  of  the  approximation  chart.  Note  that  every  point  in  the 
approximation  chart  corresponds  to  a  specific  distribution.  That  point  closest  to  the 
sample  data  locus  end  point  is  chosen  as  the  best  approximation  to  the  PDF  underlying 
the  random  data.  This  closest  point  is  determined  by  projecting  the  sample  locus  end 
point  to  all  points  on  the  approximation  chart  and  selecting  that  point  whose 
perpendicular  distance  from  the  sample  point  is  the  smallest.  Once  the  PDF  underlying 
the  sample  data  is  selected,  the  shape,  location  and  scale  parameters  are  then 
approximated. 

3.3  -  Outliers 

Outliers  that  may  exist  in  a  set  of  reference  pixels  would  seriously  change  the 
statistical  distribution  of  the  set  of  data  under  examination.  Outliers  can  cause  a  problem 
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Figure  3.2  -  Approximation  Chart,  N=NorinaI.  U=Uniform,  E=Negative 
Exponential,  A=LapIace,  S=Logistic.  C=Cauchv.  V=Extreme  Value,  T=GumbeI 
tvpe-2,  G=Gamma,  P=Pareto,  W=Weibull,  L=LognormaI,  K=K-distributed, 

B=Beta,  and  SU=Su-Johnson 


in  correctly  approximating  the  PDF  underlying  a  set  of  data  by  significantly  altering  the 
set  of  linked  vectors  generated  by  the  Ozturk  algorithm.  This  has  been  illustrated  through 
an  example  in®  where  a  set  of  NR=100  reference  data,  referred  to  as  set  A,  were  generated 
from  the  Lognormal  distribution  with  shape  parameter  0.01.  Also,  another  set,  referred  to 
as  set  B,  was  formed  which  contained  97  data  points  from  set  A  and  three  data  points 


®  M.  A.  Slamani,  “A  New  Approach  to  Radar  Detection  Based  on  the  Partitioning  and  Statistical 
Characterization  of  the  Surveillance  Volume,”  Ph.D.,  Syracuse  University,  December  1994. 
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with  very  small  values  to  constitute  the  outliers  in  the  set.  The  two  sets  were  processed  by 
the  Ozturk  algorithm  and  have  their  locus  end  points  (indicated  by  the  arrow)  reproduced 
in  the  approximation  charts  of  Figures  3.3  and  3.4.  Note  how  the  end 


point  in  Figure  3.4  for  the  set  containing  outliers  (set  B)  is  far  removed  from  the 
Lognormal  PDF  from  which  97  out  of  the  100  data  points  of  set  B  were  generated. 
Investigating  the  cause,  it  was  noted  that  the  linked  vectors  in  set  B  are  smaller  than  those 
of  set  A  causing  the  locus  end  point  for  set  B  to  fall  way  below  the  end  point  for  set  A 
and,  therefore,  outside  the  confidence  ellipses.  This  is  due  to  the  fact  that  the  amplitudes 
of  the  linked  vectors  are  proportional  to  the  magnitude  of  the  standarized  data  which 
depend  on  the  mean  and  standard  deviation  of  the  set.  The  three  outliers  do  not  affect 


29 


significantly  the  mean  of  the  set  hut  do  increase  the  variance  tremendously  causing  the 
standardized  data  and,  therefore,  the  amplitudes  of  the  linked  vectors  to  become  smaller. 
In  that  example,  the  mean  and  standard  deviation  for  set  A  are  equal  to  32.54  and  4.69, 
respectively,  while  the  mean  and  standard  deviation  for  set  B  are  equal  to  31.63  and 
18.08,  respeetively.  This  example  illustrates  what  can  happen  when  three  LP  pixels  with 
small  data  values  are  misidentified  and  associated  with  a  set  of  RPs  pixels. 


Outliers  may  be  due  to  (1)  misidentified  LP  pixels  in  the  RPs  or  vise-versa,  (2)  pixels 
having  data  values  of  very  low  probability  of  occurrence,  (3)  pixels  with  non¬ 
representative  data,  and  (4)  pixels  containing  signals  from  strong  targets.  One  way  to 
identify  outliers  within  a  region  composed  of  a  set  of  1 00  reference  pixels  is  to  compute 
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the  mean  m  and  the  standard  deviation  s  within  that  region  and  designate  as  an  outlier  any 
pixel  whose  data  value  falls  outside  the  interval  [m-ks,m+ks]  where  k  is  an  empirical 
parameter  (usually  between  1 .5  and  3)  to  be  set  by  the  user.  In  our  case  k  is  set  equal  to  2. 

3.4  -  Test  pixel  selection 

In  order  to  approximate  the  PDF  of  a  given  RPs  patch,  a  set  of  test  pixels  is 
chosen  in  that  RPs  patch.  Also,  100  reference  pixels  are  identified  for  each  test  pixel  that 
belong  to  the  selected  RPs  patch  and  are  the  closest  to  the  test  pixel  in  order  to  obtain  the 
approximating  PDF  for  each  test  pixel.  This  will  provide  information  on  how  the  data  is 
distributed  in  the  selected  RPs.  Thus  it  is  important  to  know  how  to  select  the  test  pixels 
and  their  corresponding  reference  pixels.  Test  and  reference  pixels  selection  involves 
three  steps: 

1.  A  RPs  patch  is  chosen  from  among  the  declared  RPs  patches.  This  can  be  done 
automatically  by  the  program  since  at  this  stage  every  RPs  patch  has  been  labeled 
with  a  unique  number. 

2.  A  set  of  NT  test  pixels  is  then  chosen  in  the  RPs  patch  being  processed  where  the 
value  of  NT  depends  upon  the  extent  to  which  the  patch  needs  to  be  characterized. 
Note  that  any  pixel  in  the  RPs  patch  can  be  a  test  pixel.  The  best  choice  for  the  test 
pixels  is  when  they  are  evenly  spread  throughout  the  entire  area  of  the  RPs  patch. 

3.  Finally,  for  each  test  pixel,  a  set  of  reference  pixels  is  selected.  The  reference  pixels 
must  be  in  the  same  RPs  patch  as  the  test  pixel  and  should  be  the  closest  to  it  because 
of  the  assumption  that  the  reference  pixels  are  representative  of  the  test  pixel. 

In  order  to  select  the  reference  pixels  for  a  given  test  pixel,  the  program  starts  by 
centering  a  3x3  mask  around  the  test  pixel  and  choosing  as  reference  pixels  those 
neighboring  pixels  within  the  mask  that  are  declared  to  be  in  the  same  RPs  patch  as  the 
test  pixel.  If  the  desired  number  of  reference  pixels  are  not  obtained,  the  program 
increases  the  size  of  the  mask  by  adding  one  row  and  one  column  to  each  boundary  of  the 
3x3  mask.  This  results  in  a  5x5  mask  where  only  the  pixels  in  the  augmented  rows  and 
columns  need  to  be  examined.  The  process  of  adding  one  row  and  one  column  to  each 
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boundary  of  the  previous  mask  continues  until  the  desired  number  of  reference  pixels 
have  been  obtained. 

3.5  -  Strategy  to  Sub-Patch  Investigation  Using  the  Statistical  Procedure 

The  statistical  procedure  is  applied  to  every  patch  and  subpatch  that  has  been 
declared  as  being  homogeneous  by  the  mapping  procedure.  For  each  patch  or  subpatch,  a 
set  of  test  pixels  evenly  spread  throughout  the  patch  and  their  100  closest  reference  pixels 
are  first  selected.  Let  each  set  of  100  pixels  be  referred  to  as  a  tile.  As  explained  in  the 
previous  Section,  note  that  the  sets  of  100  reference  pixels  are  chosen  to  be  disjoint,  the 
closest  to  and  belonging  to  the  same  patch  as  their  respective  test  pixels.  As  shown  in  the 
block  diagram  of  Figure  3.5,  the  statistical  procedure  consists  of  four  steps  that  are 
performed  as  follows: 

1.  Using  the  goodness  of  fit  test  of  the  Ozturk  algorithm,  a  Gaussianity  check  is 
performed  on  every  tile  to  ensure  whether  the  data  in  the  tile  are  Gaussian  or  not.  This 
results  on  every  patch  having  its  tiles  labeled  as  either  Gaussian  or  non-Gaussian. 
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Figure  3.5  -  Statistical  Procedure 
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2.  Existing  outliers  are  located  in  those  tiles  declared  as  non-Gaussian. 

3.  For  every  non-Gaussian  declared  tile,  pixels  with  outliers  are  excised  from  the  tile  and 
replaced  with  the  closest  pixels  to  the  tile  whose  data  are  not  outliers.  Next,  the 
Gaussianity  check  is  performed  once  again  as  in  step  (1).  At  the  end  of  this  step,  each 
patch  will  have  tiles  declared  as  Gaussian  or  non-Gaussian  and  pixels  with  outliers. 

4.  Using  the  Ozturk  algorithm,  the  (U,V)  coordinates  of  the  locus  end  point  is  obtained 
for  every  tile  declared  as  non-Gaussian  in  step  (3).  Next,  a  check  is  made  to  ensure 
whether  or  not  the  data  of  the  set  of  tiles  which  constitutes  a  sub-patch  can  fit  within  a 
confidence  ellipse.  This  is  done  by  first  declaring  a  subpatch  every  set  of  contiguous 
non-Gaussian  tiles.  Then,  computing  the  average  (Uav,Vav)  coordinates  of  all  test 
pixels  of  the  same  subpatch  and  getting  its  best  approximating  PDF  and  the 
corresponding  confidence  ellipse.  Finally,  a  check  is  made  whether  all  (U,V) 
coordinates  of  the  test  pixels  are  within  the  confidence  ellipse  of  the  average 
coordinates  (Uav,Vav).  If  not,  the  tiles  are  regrouped  so  that  all  (U,V)  coordinates  for 
each  group  of  tiles  can  fit  within  the  same  ellipse  as  their  corresponding  average 
(Uav,Vav)  coordinates.  Each  group  forms  then  a  subpatch. 

When  the  statistical  procedure  ends,  each  pixel  in  every  patch  is  declared  as 
Gaussian,  non-Gaussian,  or  outlier.  In  addition,  existing  subpatches  formed  by  the  sets  of 
contiguous  tiles,  whose  (U,V)  coordinates  fit  under  the  same  confidence  ellipse  as  their 
corresponding  average  (Uav,Vav)  coordinates,  are  detected. 
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4 


Expert  System  Shell  IPUS 


4.1  -  Introduction 

The  expert  system  shell  used  to  control  the  different  stages  of  the  mapping  and 
indexing  procedures  is  the  Integrated  Processing  and  Understanding  of  Signals  (IPUS) 
developed  jointly  by  the  University  of  Massachussets  and  Boston  University.  IPUS 
architecture  utilizes  the  fact  that  signal  processing  theories  supply  the  system  designer 
with  a  signal  processing  algorithm  (SPA)  that  has  adjustable  control  parameters.  Each 
instance  corresponding  to  a  particular  set  of  fixed  values  for  the  control  parameters  is 
referred  to  as  an  SPA  instance.  The  IPUS  architecture  is  designed  to  search  for 
appropriate  SPA  instances  to  be  used  in  order  to  accurately  model  the  unknown 
environment.  The  search  is  initiated  by  detection  of  a  discrepancy  at  the  output  of  a  given 
SPA  due  to  the  fact  that  the  signal  being  monitored  by  the  SPA  does  not  satisfy  the 
requirements  of  the  SPA  instance.  Once  a  discrepancy  has  been  detected,  a  diagnosis 
procedure  is  used  to  identify  the  source  of  the  distortion  that  may  have  led  to  the 
discrepancy.  Then,  either  parameters  of  the  same  SPA  can  be  readjusted  or  a  different 
SPA  can  be  chosen  to  reprocess  the  data. 

In  our  case  of  interest,  each  block  in  the  different  stages  of  the  A'SCAPE 
processor  consists  of  an  SPA  and  SPA  instances.  Rules  have  been  developed  which 
enable  the  detection  of  discrepancies  at  the  output  of  the  SPAs,  and  identification  of 
different  possible  distortion  sources  that  would  cause  the  discrepancies.  Note  that  the 
arrows  connecting  different  blocks  in  Figures  1.4,  2.8,  and  3.5  are  bi-directional.  This  is 
to  allow  for  a  system  to  assess  its  decisions,  correct  any  discrepancies,  and  adapt  to  any 
changes  in  the  environment  being  monitored. 

4.2  -  IPUS 

The  IPUS  architecture  has  evolved  from  research  on  the  design  of  a  sound 
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understanding  system'”.  The  goal  of  such  a  system  is  to  identify  the  origins  of  various 

sound  sources  (such  as  telephones,  vacuum  cleaners,  crying  infants,  etc.).  The  complexity 

of  the  sound  understanding  problem  arises  because  of  two  factors: 

1.  The  need  to  process  a  tremendous  variety  of  signal  types  due  to  the  situation 
dependent  nature  of  the  input.  For  example,  not  only  may  the  input  of  a  sound 
understanding  system  include  different  types  of  signals,  such  as  narrow-band, 
impulsive,  harmonic  signals,  but  may  also  include  various  combinations  of  these 
signals. 

2.  The  need  to  change  processing  goals  in  a  context  dependent  way.  For  example,  the 
goal  of  a  signal  understanding  system  might  be  to  respond  to  either  the  sounds  of  an 
infant  or  a  ringing  telephone  and  to  ignore  other  sound  sources.  If  an  infant  sound  is 
detected,  the  system  would  then  ignore  the  telephone  and  would  switch  its  main  goal 
to  determining  whether  the  infant  is  laughing,  crying  or  choking. 

These  two  factors  also  arise  in  the  image  processing  problem  addressed  in  this  report. 

Specifically,  complexity  is  encountered  because  of: 

1.  The  need  to  process  a  tremendous  variety  of  signal  types  due  to  the  situation 
dependent  nature  of  the  input.  For  example,  the  PDFs  of  the  random  data  in  the  pixels 
may  be  Gaussian,  Weibull,  K-distributed,  etc.,  with  various  values  for  the  scale, 
location,  and  shape  parameters. 

2.  The  need  to  change  processing  goals  in  a  context  dependent  way.  For  example,  the 
usual  operational  mode  in  imaging  is  the  enhancement  and  segmentation  of  regions.  If 
a  region  has  suspicious  characteristics  (e.g.,  isolated  tiny  region  that  may  be  a  target, 
or  signs  of  a  tumor  in  medical  imaging)  the  operational  mode  should  change  into 
isolating  the  area  and  investigating  its  nature. 


H.  Nawab,  V.  Lesser,  et  al.,  Integrated  Processing  and  Understanding  of  Signals  in  Symbolic  and 
Knowledge-Based  Signal  Processing.  Edited  by  A.V.  Oppenheim  and  H.  Nawab,  Prentice  Hal, 
Ennglewood  RPsiffs,  N.J.,  1991. 
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The  algorithms  employed  in  the  IPUS  architecture  to  identify  the  mathematical 
models  for  all  operations  of  A’ SCAPE  (e.g.,  partitioning,  detecting  edges,  approximating 
PDF)  are  referred  to  as  signal  processing  algorithms  (SPA's). 

Because  of  the  two  factors  mentioned  previously,  it  is  very  difficult,  or  even 
impossible,  to  design  a  single  mathematically  derived  signal  processing  algorithm  that 
can  be  applied  to  all  possible  input  signals  to  produee  the  desired  information  for  each 
input.  To  deal  with  such  complexities,  the  approach  taken  in  the  IPUS  architecture  is  for 
the  signal  understanding  system  to  have  access  to  a  "data-base"  of  mathematically 
derived  algorithms.  For  the  imaging  problem,  examples  are  partitioning  algorithm,  edge 
enhancement  and  detection  algorithms,  goodness  of  fit  test  algorithm,  and  the  PDF 
approximation  algorithm.  This  data  base  is  indexed  by  the  type  of  assumptions  made 
about  the  input  signal  and  the  type  of  output  information  desired  in  accordance  with  the 
eurrent  goals  of  the  signal  understanding  system. 

For  example,  it  may  be  assumed  that  the  region  under  investigation  is  Gaussian.  A 
goodness  of  fit  test  algorithm  would  be  applied  to  determine  whether  the  data  is 
statistically  consistent  with  the  Gaussian  assumption.  If  the  Gaussian  assumption  is  not 
rejected,  then  the  desired  output  information  would  be  the  location  and  scale  parameters. 

The  IPUS  architecture  utilizes  the  fact  that  signal  processing  theories  often  supply 
a  system  designer  with  a  signal  processing  algorithm  that  has  adjustable  control 
parameters  (sampling  period  of  data  samples,  number  and  location  of  reference  pixels, 
etc.).  SPA  denotes  a  data  base  of  SPA  "instances",  each  instance  corresponding  to  a 
particular  set  of  fixed  values  for  the  control  parameters.  The  IPUS  architecture  is 
designed  to  seareh  for  appropriate  SPA  instances  to  be  utilized  in  partieular  situations  in 
order  to  accurately  model  the  unknown  environment. 

Two  basic  approaches  for  carrying  out  the  signal  processing  are: 

1.  Process  the  incoming  signal  with  all  the  SPA-instances  that  are  potentially  relevant  to 
the  entire  class  of  possible  input  signals  in  the  application  domain  and  then  choose  the 
output  data  that  has  the  most  eonsistent  interpretation.  This  approach  requires  vast 
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amounts  of  signal-processing  output  data  to  be  examined  by  the  higher  level 
interpretation  processes. 

2.  Proeess  the  incoming  signal  with  one  or  a  small  number  of  the  possibly  relevant  SPA 
instances,  then  use  some  mechanism  to  recognize  if  incorrect  processing  has  taken 
place.  This  is  followed  by  determining  the  nature  of  the  incorrect  processing  through 
a  diagnostic  reasoning  process,  and  finally  changing  the  parameter  settings  of  the 
SPA  with  the  goal  of  obtaining  an  SPA-instance  which  is  appropriate  for  the 
processing  of  the  input  signal.  The  SPA-instance  with  adjusted  control  parameter 
settings  is  then  used  to  reprocess  the  input  signal. 

In  order  to  select  appropriate  values  for  the  SPA  control  parameters,  the  system 
must  consider  the  current  system  goals  as  well  as  knowledge  about  certain  characteristics 
of  the  particular  input  signal.  This  leads  to  the  dilemma  that  choosing  the  appropriate 
control  parameter  values  requires  knowledge  about  the  signal,  but  this  knowledge  can 
only  be  obtained  by  first  processing  the  signal  with  an  algorithm  with  appropriate  control 
parameter  setting.  The  IPUS  arehitecture  uses  an  iterative  technique  for  converging  to  the 
appropriate  control  parameter  values.  The  technique  begins  by  using  the  best  available 
guess  for  the  SPA  control  parameters  values.  The  SPA  instance  output  is  then  analyzed 
through  a  discrepancy  detection  mechanism  for  indicating  the  presence  of  distorted  SPA 
output  data.  A  diagnosis  is  then  performed  for  mapping  the  detected  discrepancies  to 
distortion  hypotheses.  A  signal  reprocessing  phase  then  proposes  a  new  set  of  values  for 
the  control  parameters  of  the  SPA  with  the  goal  of  eliminating  the  hypothesized 
distortions.  The  SPA  instance  corresponding  to  the  new  control  parameter  values  is  then 
used  to  reprocess  the  input  signal.  The  output  from  the  reproeessing  once  again  undergoes 
discrepancy  detection  and  if  necessary  is  followed  by  diagnosis,  signal  reprocessing 
planning,  and  further  reprocessing  of  the  input  signal. 

The  signal  data  and  the  interpretation  hypotheses  derived  from  that  data  are  stored 
on  a  blackboard  with  hierarchically  organized  information  levels.  The  hypotheses  on  the 
blackboard  fall  into  two  basic  categories:  hypotheses  posted  to  explain  the  signal  data  and 
hypotheses  posted  to  specify  expectations  about  the  nature  of  the  signal  data.  The 
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inferencing  on  the  blackboard  is  performed  by  different  knowledge  sources  (KS's)  for 
tasks  such  as  discrepancy  detection,  diagnosis  and  reprocessing,  and  data  interpretation. 
These  tasks  are  presented  in  the  following  sections. 

4.2.1  -  Discrepancy  Detection 

Ideally,  application  of  an  SPA  instance  to  input  data  results  in  undistorted  output 
data.  However,  if  the  control  parameters  of  the  SPA  instance  are  not  appropriately 
chosen,  distorted  output  data  may  result.  The  key  to  discrepancy  detection  is  the  ability  to 
recognize  and  Classify  discrepancies  due  to  distortion  introduced  by  the  SPA  instance. 
Three  types  of  discrepancies  are  possible: 

1  -  The  first  type  of  discrepancy  is  referred  to  as  a  violation.  A  violation  occurs 
when  the  SPA  output  data  implies  the  presence  of  a  signal  that  is  not  a  member  of  the 
allowable  Class  of  input  signals.  For  example,  disturbances  arising  from  pixels  in  a  clear 
region  are  always  modeled  as  Gaussian  processes  because  of  the  expectation  that 
background  noise  is  Gaussian.  Suppose  that  the  output  data  from  an  SPA  instance  implies 
that  the  disturbance  from  a  pixel  in  the  Clear  region  is  non-Gaussian.  This  constitutes  a 
violation  type  of  discrepancy. 

2  -  The  second  type  of  discrepancy  is  referred  to  as  a  conflict.  A  conflict  occurs 
when  the  current  SPA  output  data  is  inconsistent  with  expectations  arising  from 
interpretations  of  past  data.  There  are  two  types  of  conflicts  depending  upon  whether  all, 
or  only  a  portion,  of  the  current  SPA  output  data  is  inconsistent.  For  an  example  of  the 
first  type  of  conflict,  suppose  previous  SPA  output  data  arose  from  disturbances  in  the 
clear  region  while  current  SPA  output  data  is  arising  from  disturbances  in  a  patch.  A 
conflict  of  the  first  type  occurs  if  all  of  the  current  SPA  output  data,  such  as  an  increase  in 
variance  and  non-Gaussianity  of  the  data,  conflict  with  previous  interpretations  from  the 
clear  region;  For  an  example  of  the  second  type  of  conflict,  suppose  that  previous  SPA 
output  data  has  resulted  in  the  interpretation  that  the  disturbance  is  from  the  Clear  region. 
This  might  be  implied  by  the  SPA  output  data  indicating  Gaussian  statistics,  zero  mean, 
and  a  variance  level  in  the  range  of  the  background  noise.  A  conflict  of  the  second  type 
occurs  when,  even  though  Gaussian  statistics  are  confirmed  by  the  current  SPA  output 
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data,  they  also  reveal  that  the  mean  is  no  longer  zero  and  the  variance  level  has  increased 
significantly.  This  could  happen  if  the  disturbance  is  now  coming  from  a  Gaussian  patch 
where  the  data  is  highly  correlated. 

3  -  The  third  type  of  discrepancy  is  referred  to  as  a  fault.  A  fault  can  also  arise  in 
two  different  ways.  The  first  kind  occurs  when  two  or  more  different  SPA's  that  are 
applied  to  the  same  data  result  in  different  output  interpretations.  The  other  kind  occurs 
when  two  or  more  instances  of  a  single  SPA  (i.e.,  the  same  SPA  with  different  parameter 
values)  result  in  different  interpretations  when  applied  to  the  same  data.  An  example  of 
the  first  kind  of  fault  would  be  the  situation  where  SPA  #  1,  a  power  level  detector, 
indicates  a  power  level  consistent  with  the  background  noise  while  SPA  #  2,  Ozturk's 
distribution  identification  algorithm,  indicates  a  non-Gaussian  distortion.  This  is  a  fault 
because  the  background  noise  is  assumed  to  be  Gaussian.  An  example  of  the  second  kind 
of  fault  would  be  the  situation  where  use  of  Ozturk's  algorithm  based  on  1 00  and  200 
samples  from  the  same  patch  result  in  a  different  interpretation. 

4.2.2  -  Diagnosis  and  Reprocessing 

When  the  signal  being  monitored  does  not  satisfy  the  requirements  of  the  SPA 
instance,  the  output  of  the  SPA  is  distorted  resulting  in  a  discrepancy.  Once  a  discrepancy 
has  been  detected,  a  diagnosis  procedure  is  used  to  identify  the  distortion  that  may  have 
led  to  the  discrepancy.  Knowing  the  distortion,  then  either  the  appropriate  parameters  of 
the  same  SPA  can  be  adjusted  or  a  different  SPA  chosen  to  reprocess  the  data.  In  a  sense, 
the  diagnosis  procedure  maps  symptoms  (discrepancies)  to  hypothesized  underlying 
causes  (distortions).  For  example,  assume  the  sample  mean  of  a  patch  is  repeatedly  being 
evaluated  by  processing  100  samples  at  a  time.  Although,  the  first  eight  trials  result  in 
values  close  to  zero,  the  nineth  trial  produces  a  large  negative  value  for  the  mean.  This 
represents  a  conflict  of  the  first  kind.  The  diagnosis  procedure  may  surmise  that  the 
conflict  may  be  due  to  the  presence  of  one  or  more  outliers.  Consequently,  the 
reprocessing  procedure  concludes  that  the  data  from  the  ninth  trial  should  be  reprocessed 
using  a  median  detector. 
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4.2.3  -  Interpretation  Process 


The  interpretation  process  is  a  search  through  a  space  of  sets  of  interpretation 
models  for  modeling  signals.  When  a  possible  combinatorial  explosion  in  interpretation 
models  does  not  exist,  the  interpretation  process  may  be  viewed  as  just  a  straight  forward 
classification  process.  Otherwise,  the  search  must  be  carried  out  as  a  constructive 
problem  solving  process.  The  IPUS  architecture  employs  the  constructive  problem 
solving  approach  which  reduces  to  the  classification  approach  in  the  absence  of  a 
combinatorial  explosion. 

Constructive  problem  solving  techniques  must  be  used  when  the  set  of  possible 
solutions  is  too  large  to  be  enumerated.  For  example,  although  the  set  of  PDF  types  is 
finite  in  the  radar  problem,  there  are  an  infinity  of  different  PDF's  possible  because  of  the 
infinity  of  values  that  can  be  assumed  by  the  scale,  location,  and  shape  parameter. 
Consequently,  constructive  problem  solving  is  needed  to  approximate  the  underlying 
probability  distribution  of  the  data. 

4.3  -  Use  of  IPUS  in  A’SCAPE 

The  IPUS  architecture  is  suitable  when  a  single  SPA-instance  cannot  correctly 
process  all  the  input  signals  that  can  potentially  arise  in  a  signal  understanding 
application.  In  the  scene  understanding  problem,  the  variety  of  average  magnitude  levels 
in  patches  and  probability  distributions  underlying  the  data  along  with  the  different  tasks 
to  be  carried  out  in  A’SCAPE  during  the  monitoring  the  environment  (mapping, 
statistical  processing,  and  pixel  indexing)  necessitates  more  than  one  SPA-instance. 
Hence,  IPUS  is  suitable  for  the  scene  understanding  problem.  Next,  use  of  IPUS  to 
control  the  mapping  and  statistical  procedures  is  discussed. 

4.4  -  IPUS  and  the  Mapping  procedure 

Observations  are  made  next  on  the  different  control  parameters  to  set  the  stage  for 
the  design  of  the  different  rules  to  enable  IPUS  to  control  the  mapping  procedure. 
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4.4.1  -  Observations  on  Setting  the  Control  Parameters 


In  this  section,  different  effects  of  the  control  parameters  are  discussed.  Note  first 
that  the  intervals  for  the  allowable  values  of  the  control  parameters  are  given  in  Chapter 
2.0  and  are  equal  to, 

0%  <  LPQP  <  100% 

5<NCQ<8  (4.1) 

1  <  NCC  <  4 

Recall  that  LPQP  represents  the  fraction  of  LP  pixels  in  the  quantized  volume.  It 
is  used  to  determine  the  threshold  q  for  which  all  pixels  with  data  amplitudes  below  q  are 
identified  as  LP  and  all  pixels  with  data  amplitudes  above  q  are  identified  as  RPs  in  the 
quantized  volume. 

Also,  NCQ  is  the  minimum  number  of  neighboring  pixels  in  the  quantized 
volume  required  to  be  identified  as  RPs  pixels  for  a  test  pixel  to  be  declared  as  a  RPs 
pixel  in  the  first-corrected  volume.  Finally,  NCC  is  the  minimum  number  of  neighboring 
pixels  in  the  first-corrected  volume  required  to  be  identified  as  RPs  pixels  for  a  test  pixel 
to  be  declared  as  a  RPs  pixel  in  the  second-corrected  volume. 

LPCQP  and  LPCCP  are  computed  parameters  which  represent  the  LP  percentages 
in  the  first  and  second-corrected  volumes,  CQV  and  CCV,  respectively. 

Define  LPQPt  to  be  the  true  value  for  the  fraction  of  LP  pixels  in  the  generated 

scene. 

As  explained  in  Chapter  2  the  mapping  processor  begins  by  setting  a  threshold 
that  results  in  a  specified  fraction  of  LP  pixels  equal  to  LPQP.  The  mapping  processor 
iterates  until  it  is  satisfied  that  the  latest  scene  is  consistent  with  the  last  specified  value  of 
LPQP.  When  the  iteration  process  ends,  it  is  assumed  that 

LPQP  s  LPQPt  (4.2  ) 
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4.4.1. a  -  Observations  on  the  Setting  of  LPOP 

1.  Setting  LPQP  much  smaller  than  LPQPt:  Many  pixels  have  data  amplitudes  larger 
than  the  threshold  resulting  in  a  large  number  of  LP  pixels  being  declared  as  RPs 
pixels  in  the  quantized  volume. 

2.  Setting  LPQP  much  larger  than  LPQPt:  RPs  patch  pixels  may  be  misclassified  due 
to  the  fact  that  some  RPs  patches  have  data  amplitude  values  below  the  threshold. 
This  results  in  many  RPs  pixels  being  identified  as  LP  pixels  in  the  quantized  volume. 

Conclusion:  The  threshold  is  always  set  very  low  at  the  begirming  so  that  LP 
information  is  gained  as  the  process  iterates.  Because  (1)  the  objective  of  the  mapping 
procedure  is  to  separate  between  LP  and  RPs  patches,  (2)  the  average  power  of  the  LP  is 
the  lowest  among  all  regions,  and  (3)  the  threshold  is  set  adaptively  by  the  assessment 
stage,  the  threshold,  controlled  by  the  assessment  stage,  is  raised  adaptively  until 
LPQP  s  LPQPt. 

4.4.1. b  -  Observations  on  the  Setting  of  NCO 

Recall  that  NCQ  controls  which  test  pixels  in  the  first-corrected  volume  are  to  be 
declared  as  RPs.  NCQ  is  said  to  be  large  when  its  value  approaches  8  and  small  when  its 
value  approaches  5.  The  following  observations  relative  to  NCQ  take  into  consideration 
that  the  initial  setting  of  LPQP  is  low  and  then  is  increased  until  LPQP  approximates  the 
true  value  LPQPt.  Depending  on  the  setting  of  LPQP  with  respect  to  LPQPt,  four  cases 
exist: 

1.  LPQP  much  smaller  than  LPQPt: 

a)  Setting  NCQ  small:  In  this  case,  because  many  RPs  declared  pixels  exist  in 
the  quantized  volume  due  to  the  low  threshold,  small  NCQ  results  in  the 
building  of  a  multitude  of  RPs  patches  which  are  likely  to  be  so  close  that  they 
form  a  single  big  RPs  patch  in  the  first-corrected  volume. 

b)  Setting  NCQ  high:  Here,  even  though  many  RPs  declared  pixels  exist  in  the 
quantized  volume  due  to  the  low  threshold,  high  NCQ  results  in  the  building 
of  fewer  RPs  patches  than  when  NCQ  is  small.  This  is  due  to  the  fact  that 
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there  must  be  at  least  NCQ  RPs  pixels  neighboring  the  test  pixel  in  the 
quantized  volume,  where  NCQ  is  large,  in  order  for  the  test  pixel  to  be 
declared  as  a  RPs  pixel  in  the  first-corrected  volume.  In  this  case,  corrections 
are  made  and  some  of  the  pixels  previously  declared  as  RPs  pixels  in  the 
quantized  volume  are  now  declared  as  LP  pixels  in  the  first-corrected  volume. 

2.  LPQP  close  to  LPQPt:  When  LPQP  is  close  to  its  true  value,  the  threshold  is  high 
enough  to  separate  between  the  LP  region  and  RPs  patches.  With  either  small  or  large 
values  for  NCQ,  the  RPs  regions  are  well  approximated.  In  this  case,  the  choice  of 
NCQ  affects  the  classification  of  the  inner  pixels  of  the  RPs  regions.  This  results 
because,  even  though  the  data  amplitudes  of  RPs  pixels  are  higher  than  those  of  the 
LP  pixels,  in  general,  some  RPs  pixels  with  data  amplitudes  lower  than  those  of  the 
highest  LP  data  values  exist  and  may  be  lower  than  the  threshold. 

a)  Setting  NCQ  small:  All  test  pixels  in  the  quantized  volume  that  have  at  least 
NCQ  neighboring  pixels  are  declared  as  RPs  pixels  in  the  first-corrected 
volume.  Small  NCQ  helps  to  correctly  classify  the  inner  RPs  pixels.  However, 
note  that  small  NCQ  also  results  in  misclassifying  LP  pixels  that  are 
surrounded  by  at  least  NCQ  declared  RPs  pixels. 

b)  Setting  NCQ  high:  Every  test  pixel  must  have  a  large  number  of  neighboring 
RPs  declared  pixels  in  the  quantized  volume  for  it  to  be  declared  as  a  RPs 
pixel  in  the  first-corrected  volume.  This  causes  the  procedure  to  misclassify 
some  of  the  inner  RPs  pixels  when  too  many  of  the  neighboring  pixels  have 
their  data  amplitudes  falling  below  the  threshold.  In  this  case,  the  identified 
RPs  regions  are  not  homogeneous  and  contain  LP  declared  "holes". 

Conclusion:  The  value  of  NCQ  should  be  chosen  as  large  as  possible  at  the 
beginning  of  the  iterative  process  when  the  threshold  is  set  very  low  to  correctly 
reclassify  the  maximum  number  of  LP  pixels  misidentified  at  the 
thresholding/quantization  stage.  When  the  threshold  reaches  a  level  where  it  is  close  to  its 
convergence  value,  NCQ  should  then  be  chosen  small  to  avoid  non-homogeneous  RPs 
regions. 
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4.4.1.C  -  Observations  on  the  Setting  of  NCC 


Because  NCQ  truncates  the  boundaries  of  the  RPs  regions,  NCC  is  used  to 
augment  the  edges  of  the  RPs  declared  regions.  NCC  is  said  to  be  large  when  its  value 
approaches  4  and  small  when  its  value  approaches  1.  In  the  following  discussion  it  is 
assumed  that  the  conclusions  previously  reached  on  the  settings  of  LPQP  and  NCQ  are 
taken  into  consideration  so  that  LPQP  is  initially  set  low  to  be  increased  until  it 
approaches  its  true  value  LPQPt,  while  NCQ  is  initially  set  to  a  large  value,  to  be 
decreased  as  LPQP  approaches  its  true  value.  Four  cases  are  then  identified: 

1.  LPQP  much  smaller  than  LPQPt  and  NCQ  large:  Because  NCQ  is  set  large,  many 
RPs  edge  pixels  are  misclassified  and  associated  with  the  LP  region. 

a)  Setting  NCC  small:  When  NCC  is  set  small,  many  of  the  edge  pixels  are 
correctly  reclassified  from  LP  pixels  to  RPs  pixels  in  the  second-corrected 
volume. 

b)  Setting  NCC  large:  In  this  case,  only  a  few  misclassified  RPs  edge  pixels  are 
correctly  reclassified  in  the  second-corrected  volume. 

2.  LPQP  close  to  LPQPt  and  NCQ  small:  Because  NCQ  is  small,  only  a  few  RPs  edge 
pixels  are  associated  with  the  LP. 

a)  Setting  NCC  small:  Small  NCC  causes  not  only  RPs  edge  pixels  to  be 
recovered  but  also  LP  pixels  to  be  misclassified  in  the  second-corrected 
volume. 

b)  Setting  NCC  large:  In  this  case,  most  of  the  RPs  edge  pixels  are  correctly 
classified  in  the  second-corrected  volume  and  only  few  LP  pixels  are 
misclassified  as  RPs  pixels. 

Conclusion:  NCC  results  in  the  recovery  of  RPs  edge  pixels  and  the 
misclassification  of  some  LP  pixels  close  to  the  RPs  edge  pixels.  In  order  to  maximize 
recovery  of  the  RPs  edge  pixels  and  minimize  the  misclassification  of  LP  pixels,  NCC 
should  be  set  small  when  NCQ  is  set  large  in  order  to  recover  many  of  the  RPs  edge 
pixels  that  were  lost  in  the  first-correction.  On  the  other  hand,  NCC  should  be  set  large 
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when  NCQ  is  set  small  because,  in  this  case,  only  a  few  RPs  edge  pixels  need  to  be 
recovered. 

4.4.2  -  Resolution  of  Discrepancies 

In  this  section,  rules  are  developed  to  enable  the  resolution  of  discrepancies.  The 
assessment  stage  of  the  mapping  procedure  consists  of  comparing  at  each  step  of  the 
iteration  the  value  for  LPCCP  with  the  corresponding  LPQP.  When  LPCCP  is  not 
sufficiently  close  to  LPQP,  the  assessment  stage  is  said  to  fail.  This  initiates  the 
discrepancy  detection  stage.  Diagnosis  identifies  the  distortion  that  may  have  caused  the 
discrepancy  and  adjusts  one  or  more  of  the  mapping  control  parameters  for  reprocessing 
of  the  data. 

The  strategy  behind  the  iterative  process  of  the  mapping  procedure  employs  two 
stages.  In  the  first  stage,  referred  to  as  the  threshold  approximation  stage,  LPQP  is  varied 
iteratively  by  the  mapping  processor  until,  as  explained  later,  it  is  expected  that  LPQP  is 
within  10%  of  its  true  value  LPQPt.  The  second  stage,  referred  to  as  the  threshold  fine- 
tuning  stage,  consists  of  iteratively  varying  LPQP  until  it  converges  to  within  1%  of  the 
last  computed  value  for  LPCCP.  The  two  stages  are  now  discussed  in  detail. 

4.4.2.a  -  Discrepancies  in  the  Threshold  approximation  stage 

During  this  stage,  two  sets  of  SPA  instances  are  used  on  the  same  data  of  the 
surveillance  volume.  For  both  sets  LPQP  and  NCC  are  the  same  whereas  NCQ  is  equal  to 
7  for  one  set  and  8  for  the  other. 

In  order  to  understand  how  the  rules  are  set  in  this  stage,  it  is  important  to 
subdivide  the  problem  into  two  cases.  Case  I  describes  the  situation  where  the  LP  and 
RPs  patches  have  histograms  with  large  area  of  overlapping  whereas  case  II  describes  the 
situation  where  the  LP  and  RPs  patches  have  histograms  with  small  area  of  overlapping 
to  non-overlapping. 

Case  I  -  LP  and  RPs  whose  histograms  have  large  area  of  overlap 

Recall  that  NCQ  is  used  to  recognize  the  RPs  patches  in  the  surveillance  volume. 
First  consider  the  situation  where  LPQP  approximately  equals  LPQPt.  Here  the  threshold 
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is  such  that  it  is  possible  to  do  a  good  job  of  separating  between  the  LP  region  and  RPs 
patches.  Steps  are  then  taken  to  correct  misclassified  LP  and  RPs  data.  Note  that 
misclassifications  are  due  to  large  LP  data  exceeding  the  threshold  and  small  RPs  data 
falling  below  the  threshold.  At  this  point,  setting  NCQ  to  7  and  8,  respectively,  results  in 
very  close  values  for  LPCQP  and  LPCCP  due  to  the  facts  that  (1)  the  two  masks  are  very 
similar  (NCQ=8  requires  that  8  neighboring  pixels  be  declared  RPs  in  the  quantized 
volume  for  a  test  pixel  in  the  first-corrected  volume  to  be  declared  RPs  whereas  NCQ=7 
requires  that  7  neighboring  pixels  be  declared  RPs  in  the  quantized  volume  for  a  test  pixel 
in  the  first-corrected  volume  to  be  declared  RPs),  (2)  only  a  few  pixels  are  misclassified 
in  the  quantized  volume. 

Now  consider  that  LPQP  is  significantly  smaller  than  LPQPt.  In  this  case  many 
LP  pixels  are  misclassified  after  quantization.  Even  though  masks  with  NCQ  equal  to  7 
and  8  are  similar,  they  result  in  LPCQP  and  LPCCP  being  considerably  different  due  to 
the  fact  that  the  large  number  of  misclassified  LP  pixels  are  so  many  that  they  tend  to 
group  together.  Consequently,  changing  NCQ  from  8  to  7  simply  results  in  additional  LP 
pixels  grouping  together  to  form  additional  RPs  regions  and  more  edges.  Because  of  this, 

LPCQPUcq=7<LPCQP|hcq=8 

and  (4.3 ) 

LPCCP|f^(;;Q=7  <  LPCCP|ncQ=8  ■ 

Because  the  second-corrected  volume  represents  the  scene  where  RPs  patches  and 
their  edges  are  assumed  to  be  properly  recovered  had  the  threshold  LPQP  been  chosen 
properly,  LPCCP  tries  to  converge  to  LPQPt.  Thus,  it  is  logical  to  begin  each  iteration  by 
assigning  to  LPQP  the  latest  computed  value  of  LPCCP|ncq=8-  The  very  first  value 
assigned  to  LPQP  is  simply  a  guess.  This  value  should  be  such  that  the  threshold  is  low. 
In  all  of  our  examples,  the  first  value  of  LPQP  is  chosen  equal  to  10%. 

Using  different  scenes  with  different  values  for  LPQPt,  it  has  been  determined 
near  convergence  that  whenever  the  difference  between  LPCCP|ncq=7  LPCCP|ncq=8  is 
within  10%,  then  LPQP  is  likely  to  be  within  10%  of  its  true  value.  This  is  confirmed  in 
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Figure  4.1  where  plots  of  the  quantities  LPCCP1ncq=8  -  LPCCP1ncq=7  versus  LPQP-LPQPt 
are  shown  for  different  values  of  LPQPt.  Note  that 

-  i  -  for  different  values  of  LPQPt,  all  plots  are  such  that  near  convergence  the 
difference  LPCCP|ncq=8  -  LPCCP|ncq=7  approaches  zero  when  LPQP-LPQPt  also 
approaches  zero, 

-  ii  -  for  different  values  of  LPQPt,  when  LPQP-LPQPt  >  -10%,  the  difference 
LPCCP|ncq==8  -  LPCCP|ncq=7  <  10%. 

As  a  result  of  the  above,  the  threshold  approximation  stage  iterates  until  it  is 
satisfied  that 

LPCCP|ncq=8  -  LPCCP|ncq=7  <  1 0%.  (4.4  ) 

Case  II  -  LP  and  RPs  whose  histograms  either  have  a  small  area  of  overlap  or  do  not 
overlap  at  all 

Equation  4.4  is  mostly  valid  when  the  histograms  of  the  different  patches  overlap 
noticeably.  When  the  histograms  do  not  overlap  or  have  a  very  small  region  of  overlap,  it 
is  noticed  that  the  difference  in  Equation  4.4  is  even  less  than  10%  for  all  settings  of 
LPQP,  when  LPQP  is  significantly  smaller  or  larger  than  LPQPt  many  LP  pixels  are 
misclassified  after  quantization.  Even  though  masks  with  NCQ  equal  to  7  and  8  are 
similar,  they  result  in  LPCQP  and  LPCCP  being  considerably  different  due  to  the  fact 
that  the  large  number  of  misclassified  LP  pixels  are  so  many  that  they  tend  to  group 
together.  Consequently;  changing  NCQ  from  8  to  7  simply  results  in  additional  LP  pixels 
grouping  together  to  form  additional  RPs  regions  and  more  edges  causing  the  difference 
of  LPCCP1ncq=8  -  LPCCP|ncq=7  to  be  large.  On  the  other  hand,  when  LPQP  is  very  close  to 
LPQPt  the  threshold  is  set  so  that  it  is  possible  to  separate  between  the  LP  and  RPs,  the 
setting  of  NCQ  to  7  and  8,  respectively,  results  in  the  difference  of  LPCCP1ncq=8- 
LPCCP1ncq=7  being  very  small  due  to  the  fact  that  only  a  few  pixels  are  misclassified  in 
the  quantized  volume. 
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Using  the  observations  above,  it  is  concluded  that  in  order  to  find  the  threshold  to 
separate  between  LP  and  RPs  when  their  histograms  do  not  overlap  it  is  sufficient  to  look 
for  the  minimum  in  the  difference  LPCCP|ncq=8  -  LPCCP|ncq=7  as  LPQP  is  varied. 

In  summary,  a  guess  for  the  initial  value  of  LPQP  is  followed  by  the  execution  of 
the  mapping  procedure  using  two  different  SPA  instances.  The  outputs  of  the  two  SPA 
instances  are  compared  by  means  of  the  computed  values  of  LPCCP.  If  LPCCP|ncq=8- 
LPCCP|ncq=7  is  more  than  10%,  a  discrepancy  is  detected  and  it  is  concluded  that  the 
value  of  LPQP  differs  from  its  true  value  by  more  than  10%.  LPQP  is  then  increased  to 
the  latest  computed  value  of  LPCQP|ncq=8-  LPQP  is  varied  from  one  iteration  to  the  next 
while  NCC  is  kept  equal  to  1 .  This  choice  for  NCC  agrees  with  the  observations  made 


Figure  4.1  -  Plot  of  LPCCPI^cq,.  -  LPCCP^co^?  Versus  LPOP-LPOPt  for  Different 

Values  of  LPOPt 


previously  where  it  was  concluded  that  NCC  should  be  set  small  when  NCQ  is  set  large. 
In  this  case  NCQ  has  a  large  value  equal  to  either  7  or  8. 
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The  discrepancies  that  may  arise  in  the  threshold  approximation  stage  are  due  to 
the  fact  that  two  instances  of  the  same  SPA  result  in  different  interpretations  when 
applied  to  the  same  data.  As  defined  in  Section  4.2.1,  such  a  discrepancy  is  typified  as  a 
fault. 

Four  fault-type  discrepancies  are  readily  identified  in  the  threshold  approximation 
stage.  These  are  as  follows: 

1.  LPCCP|ncq=8  -  LPCCPUcq=7  >  10%  :  As  discussed  above,  the  goal  of  the 
threshold  approximation  stage  is  to  obtain  a  threshold  LPQP  that  is  within 
10%  of  its  true  value.  As  shown  in  Figure  4.1,  this  is  likely  only  when 
LPCCP|ncq=8  -  LPCCP1ncq=7  <  10%-  When  the  difference  between  the 
computed  thresholds  LPCCP|ncq=8  -  LPCCP|ncq=7  is  more  than  10%,  a  fault 
type  of  discrepancy  is  detected  during  the  assessment  stage.  The  diagnosis 
process  identifies  the  fact  that  LPQP-LPQPt  <  -10%  as  the  source  for  the 
distortion  causing  the  discrepancy.  The  remedy,  in  this  case,  is  to  increase  the 
value  of  LPQP  during  the  reprocessing  stage  to  the  latest  computed  value  for 
LPCCPUcq=8. 

2.  initial  LPQP  set  too  low:  In  some  cases,  when  the  initial  guess  for  LPQP  is 
too  small,  the  number  of  LP  pixels  with  data  exceeding  the  threshold  is  so 
large  that  when  corrections  are  made,  the  second-corrected  volume  results  in 
many  RPs  declared  patches  or,  in  the  worst  case,  a  single  big  RPs  patch.  This 
results  in  the  values  of  either  LPCCP|ncq=8  or  LPCCP|ncq=7  being  even  smaller 
than  LPQP.  In  this  case,  it  is  possible  to  obtain  a  value  for  the  difference 
LPCCP|ncq=8  -  LPCCP|ncq=7  that  is  smaller  than  10%.  The  IPUS  control  must 
be  suspicious  of  such  a  case  and  declare  LPQP-set-too-low  as  the  source  for 
the  distortion  causing  the  discrepancy.  The  remedy,  in  this  case,  is  to  increase 
the  initial  value  of  LPQP  during  the  reprocessing  stage.  In  our  examples,  we 
choose  to  increase  LPQP  by  10%  every  time  an  initial-LPQP-set-too-low  fault 
is  obtained. 
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3.  LPCCP|[v,cq=8  -  LPCCP|ncq=7  <  10%  even  when  LPQP  is  increased:  The 
IPUS  control  must  be  suspicious  of  such  a  case  and  declare  the  fact  that  the 
histograms  of  the  LP  and  RPs  do  not  overlap  as  the  source  for  the  distortion 
causing  the  discrepancy.  The  remedy,  in  this  case,  increase  LPQP  by  10% 
steps  until  a  minimum  value  for  LPCCP|NCQ=g  -  LPCCP|ncq=7  is  obtained. 

4.  It  is  noted  that  during  the  threshold  approximation  stage  it  is  possible  that  the 
inequality  in  Equation  4.4  will  be  satisfied  with  PLCCP|,,jcq=7  or  s  =100%. 
This  means  that  the  subpatch  with  the  smallest  average  power  occupies  100% 
of  the  RPs  patch  area.  Consequently,  the  RPs  patch  is  homogeneous.  When 
the  inequality  in  Equation  4.4  is  met  with  PLCCP|ncq=7  =100%,  there  is  no 
need  for  the  threshold  fine-tuning  stage.  This  is  because  the  initial  value  of 
PLQP  in  the  threshold  fine-tuning  stage  would  be  equal  to  100%  and  any 
more  processing  would  also  result  in  PLCCP=100%  regardless  of  the  values 
chosen  for  NCQ  and  NCC. 


Table  4.1,  summarizes  the  discrepancies  that  may  occur  during  the  threshold 
approximation  stage. 


Discrepancy 

(Fault  Type) 

Diagnosis 

Reprocessing 

LPCCP|ncq.8  -  LPCCPUcq^, 
is  more  than  1 0% 

PLQP-PLQPt>-10% 

Assign  to  PLQP  the  latest 
value  of  PLQP|ncq=8 

LPCCPI^cq^s  -  LPCCPUcq=7 
is  less  than  10%  in  the  early 
stages  of  iteration 

Initial  value  of  PLQP  is 
too  low 

Increase  PLQP  by  10%  from 
its  initial  value 

LPCCPUcq^s  -  LPCCPIncq^, 
is  less  than  10%  in  all 
iterations 

LP  and  RPs  histograms 
do  not  overlap 

Search  for  the  PLQP  that 
results  in  a  minimum  value 
LPCCPUcq=8-LPCCPUcq=7 

Inequality  in  Eq.  (VII.4-1) 
will  be  satisfied  with 

PLCCP|NCQ=7or8=100% 

The  scene  is 

homogeneous 

Declare  that  no  LP  exist.  No 
threshold  fine-tuning  is 
needed 

Table  4.1  -  Discrepancies  in  the  Threshold  Approximation  Stage 
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4.4.2.b  -  Threshold  fine-tuning  stage 


At  the  end  of  the  threshold  approximation  stage  LPQP  is  likely  to  be  within  10% 
of  its  true  value  LPQPt.  During  the  threshold  fine-tuning  stage,  LPQP  is  varied  until 
LPCCP  is  within  1%  of  the  corresponding  value  of  LPQP. 

During  the  threshold  fine-tuning  stage,  NCQ  is  lowered  to  avoid  holes  in  the  RPs 
patches  caused  by  misclassified  RPs  pixels  whose  data  values  are  lower  than  the 
threshold.  For  the  same  reason  LPQP  is  assigned  the  latest  value  of  LPCCP|ncq=7  rather 
than  the  latest  value  of  LPCCP|ncq=8-  addition,  the  value  of  NCC  is  raised  to  avoid 
misclassification  of  LP  pixels  close  to  the  RPs  edges.  These  choices  for  NCQ  and  NCC 
agree  with  observations  mentioned  in  Section  4.4. 1 . 

The  following  observations  on  LPQP,  NCQ,  and  NCC  are  necessary  to 
understand  how  these  parameters  should  be  automatically  set  in  order  for  LPCCP  to 
converge  to  within  1%  of  LPQP. 

1  -  When  LPQP  is  increased  while  both  NCQ  and  NCC  are  kept  constant,  the 
number  of  LP  pixels  in  the  quantized  volume  is  increased  and,  therefore,  both 
LPCQP  and  LPCCP  are  likely  to  increase. 

2  -  When  NCQ  is  increased  while  both  LPQP  and  NCC  are  kept  constant,  the 
requirement  on  a  test  pixel  to  be  declared  as  a  RPs  pixel  in  the  first-corrected 
volume  becomes  more  stringent  and,  therefore,  the  number  of  RPs  pixels  in  both 
corrected  volumes  are  likely  to  decrease.  This  tends  to  increase  the  number  of  LP 
pixels  causing  both  LPCQP  and  LPCCP  to  increase. 

3  -  When  NCC  is  increased  while  both  LPQP  and  NCQ  are  kept  constant,  the 
requirement  on  a  test  pixel  to  be  declared  as  a  RPs  pixel  in  the  second-corrected 
volume  becomes  more  stringent  and,  therefore,  the  number  of  RPs  pixels  in  this 
volume  is  decreased.  Thus,  the  number  of  LP  pixels  in  the  second-corrected 
volume  increases  and,  consequently,  LPCCP  increases. 

Using  the  above  observations,  the  following  strategy  is  used  by  the  assessment 
stage  to  control  the  threshold  fine-tuning  stage. 
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1  -  Because  LPQP  is  within  10%  of  its  true  value  at  the  beginning  of  the  threshold 
fine-tuning  stage,  the  threshold  is  likely  to  be  relatively  high.  Thus,  NCQ  should 
be  set  to  its  smallest  value  of  5  while,  as  needed,  NCC  should  be  incremented 
iteratively  from  its  minimum  value  of  1  up  to  its  maximum  value  of  4. 

2  -  When  the  inequality  in  Equation  4.4  is  not  satisfied,  LPQP  should  be  increased 
in  small  steps.  Otherwise,  the  iterative  process  diverges  when  the  same  rule  from 
the  threshold  approximation  stage  is  used.  The  approach  taken  in  this  work  during 
the  threshold  fine-tuning  stage  consists  of  assigning  a  value  to  LPQP  that  is  half 
way  between  its  latest  value  and  the  latest  value  of  LPCCP,  i.e. 

LPQP  =  +  LPCCP, 

3  -  The  condition  set  forward  for  ending  the  threshold  fine-tuning  stage  is  given 
by 

|LPQP  -  LPCCPI  <  1%  (4.6) 


Two  cases  are  possible  when  the  inequality  in  Equation  4.6  is  not  satisfied:  either 
LPQP  <  LPCCP  or  LPQP  >  LPCCP. 

4  -  When  the  inequality  in  Equation  4.6  is  not  satisfied  and  LPQP  <  LPCCP,  the 
control  parameters  should  be  varied  by  the  diagnosis  procedure  such  that  LPCCP 
is  decreased.  In  this  case,  NCQ  should  be  made  smaller.  If  none  of  the  allowable 
values  for  NCQ  result  in  the  inequality  in  Equation  4.6  being  satisfied,  then  LPQP 
is  varied  according  to  Equation  4.5. 

5  -  When  the  inequality  in  Equation  4.6  is  not  satisfied  and  LPQP  >  LPCCP,  the 
control  parameters  should  be  varied  by  the  diagnosis  procedure  such  that  LPCCP 
is  increased.  In  this  case,  NCC  should  be  made  larger.  If  none  of  the  allowable 
values  for  NCC  result  in  the  inequality  in  Equation  4.6  being  satisfied,  then  LPQP 
is  varied  according  to  Equation  4.5. 
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6  -  If  the  threshold  fine-tuning  stage  results  in  PLCCP=100%  at  any  iteration,  the 
threshold  fine-tuning  stage  should  end  because,  as  in  observation  2-,  any  more 
processing  will  end  with  PLCCP=100%  regardless  of  the  values  chosen  for  NCQ 
and  NCC.  This,  in  turn,  will  make  PLQP  equal  to  100%. 

Note  that  during  the  threshold  approximation  stage  note  that  only  LPQP  was 
varied  whereas,  in  the  threshold  fine-tuning  stage,  any  of  the  parameters  LPQP,  NCQ, 
and  NCC  may  be  varied. 

In  summary,  the  threshold  fine-tuning  stage  begins  by  assigning  to  LPQP  the 
latest  value  of  LPCCP|ncq=7'  Once  quantization,  first-correction,  and  second-correction 
stages  are  completed  with  pre-selected  values  for  NCQ  and  NCC,  the  assessment  stage 
diagnoses  the  results  according  to  the  strategy  discussed  above,  and,  depending  on  the 
outcome,  decides  either  that  reprocessing  is  necessary  with  adjusted  values  for  any  of  the 
LPQP,  NCQ,  and  NCC  parameters,  or  the  threshold  fine-tuning  stage  is  completed. 

At  the  end  of  each  iteration  of  the  threshold  fine-tuning  stage  it  is  expected  that 
the  computed  value  of  LPCCP  will  be  within  1%  of  LPQP.  When  the  inequality  in 
Equation  4.6  is  not  satisfied,  a  conflict  type  of  discrepancy  is  detected  based  on  the 
inconsistency  in  the  expectation  that  LPCCP  will  be  within  1%  of  LPQP.  Table  4.2 
summarizes  the  discrepancies  that  may  occur  during  the  threshold  fine-tuning  stage. 


Discrepancy 

(Conflict  Type) 

Diagnosis 

Reprocessing 

|PLQP-PLCCP1>1%  and 

PLQP<PLCCP 

Either  NCQ  or  PLQP 
are  not  well  adjusted 

Decrease  NCQ,  otherwise 
update  LPQP 

|PLQP-PLCCP1>1%  and 

PLQP>PLCCP 

Either  NCC  or  PLQP 
are  not  well  adjusted 

Increase  NCC,  otherwise 
update  PLQP 

Table  4.2  -  Discrepancies  in  the  Threshold  Fine-Tuning  Stage 
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It  has  been  determined  through  examples  that  the  initial  setting  of  NCQ=5  is 
adequate  for  the  threshold  fine-tuning  stage  to  converge  (i.e.,  in  the  different  treated 
examples,  it  was  never  necessary  to  decrease  NCQ  from  the  value  of  5). 

4.5  -  IPUS  and  the  Statistical  Procedure 

Recall  that  the  objective  of  the  statistical  procedure  is  to  separate  between 
contiguous  non-homogeneous  patches  based  on  their  data  distributions.  A  four  step 
strategy  for  the  statistical  procedure  stage  was  presented  in  Chapter  3.  The  steps  are: 

1.  Using  the  goodness  of  fit  test  of  the  Ozturk  algorithm,  a  Gaussianity  check  is 
performed  on  every  tile  to  ensure  whether  the  data  in  the  tile  are  Gaussian  or  not.  This 
results  on  every  patch  having  its  tiles  labeled  as  either  Gaussian  or  non-Gaussian. 

2.  Existing  outliers  are  located  in  those  tiles  declared  as  non-Gaussian. 

3.  For  every  non-Gaussian  declared  tile,  pixels  with  outliers  are  excised  from  the  tile  and 
replaced  with  the  closest  pixels  to  the  tile  whose  data  are  not  outliers.  Next,  the 
Gaussianity  check  is  performed  once  again  as  in  step  (1).  At  the  end  of  this  step,  each 
patch  will  have  tiles  declared  as  Gaussian  or  non-Gaussian  and  pixels  with  outliers. 

4.  Using  the  Ozturk  algorithm,  the  (U,V)  coordinates  of  the  locus  end  point  is  obtained 
for  every  tile  declared  as  non-Gaussian  in  step  (3).  Next,  a  check  is  made  to  ensure 
whether  or  not  the  data  of  the  set  of  tiles  which  constitutes  a  sub-patch  can  fit  within  a 
confidence  ellipse.  This  is  done  by  first  declaring  a  subpatch  every  set  of  contiguous 
non-Gaussian  tiles.  Then,  computing  the  average  (Uav,Vav)  coordinates  of  all  test 
pixels  of  the  same  subpatch  and  getting  its  best  approximating  PDF  and  the 
corresponding  confidence  ellipse.  Finally,  a  check  is  made  whether  all  (U,V) 
coordinates  of  the  test  pixels  are  within  the  confidence  ellipse  of  the  average 
coordinates  (Uav,Vav).  If  not,  the  tiles  are  regrouped  so  that  all  (U,V)  coordinates  for 
each  group  of  tiles  can  fit  within  the  same  ellipse.  Each  group  then  forms  a  subpatch. 

IPUS  is  primarily  needed  in  the  fourth  step  of  the  statistical  procedure  where  a 
grouping  of  contiguous  tiles  whose  (U,V)  pairs  can  fit  under  a  single  ellipse  is  to  be 
performed  to  form  subpatches.  This  step  necessitates  rules  to  enable  the  diagnosis  process 
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to  acknowledge  that  a  diserepancy  exists  that  is  eaused  by  the  distortion  due  to  the 
eontiguous  non-Gaussian  tiles  not  being  all  homogeneous.  This  distortion  is  displayed 
when  a  set  of  contiguous  non-Gaussian  tiles  cannot  have  all  of  their  (U,V)  pairs  fit  under 
the  same  ellipse. 

Then,  reproeessing  takes  place  to  find  sets  of  homogeneous  non-Gaussian  tiles 
whose  (U,V)  pairs  can  fit  under  the  same  ellipse  and  deelare  each  set  as  a  subpatch.  This 
is  done  under  the  eonstraint  that  a  minimum  number  of  ellipses  should  be  used.  Note  that 
the  discrepancy  in  this  case  is  of  eonfliet  type  because  the  expeetation  that  a  set  of 
eontiguous  non-Gaussian  tiles  should  be  homogeneous  is  not  met.  Table  4.3  summarizes 
the  discrepancy  in  the  statistical  procedure. 


Discrepancy 

(Conflict  Type) 

Diagnosis 

Reprocessing 

A  set  of  contiguous  non- 
Gaussian  tiles  eannot  have 
all  of  their  (U,V)  pairs  fit 
under  the  same  ellipse 

Not  all  of  the 
eontiguous  non- 
Gaussian  tiles  are 
homogeneous 

Find  sets  of  homogeneous 
tiles  whose  (U,V)  pairs  ean 
fit  under  the  same  ellipse  and 
declare  each  set  as  a 
subpateh. 

Table  4.3  -  Discrepancy  in  the  Statistical  Procedure 


4,6  -  Conclusion 


The  expert  system  IPUS  is  used  to  guarantee  the  convergence  of  both  the  mapping 
and  statistieal  procedures.  Through  its  rules  for  reprocessing  IPUS  eontrols  all  feed¬ 
forward  and  feed-baek  eonnections  between  the  different  blocks  of  A’ SCAPE. 
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Application  of  A*  SCAPE  to  real  IR  data 


5.1  -  Introduction 

In  this  chapter,  real  data  of  an  IR  image  are  processed  through  A’ SCAPE.  Section 
5.2,  discusses  how  the  data  is  collected.  Sections  5.3  through  5.6  present  the  results  of 
processing  the  IR  scene  through  the  different  stages  of  A’ SCAPE.  Finally,  section  5.7 
consists  of  a  conclusion  for  real  data  processing. 

5.2  -  Data  Collection 

The  data  processed  in  this  work  consist  of  real  airborne  data  of  infrared  (IR)  type 
collected  over  lake  Michigan.  As  shown  in  Figure  5.1”  (reproduced  from  the  reference 
mentioned  in  the  footnote),  the  airborne  data  is  collected  through  an  m-7  line  scanner 
which  consists  of  an  optical  telescope  with  its  narrow  field  of  view  (ground  resolution 
element)  redirected  by  a  rotating  flat  mirror  causing  the  system  to  scan  in  a  plane 
perpendicular  to  the  longitudinal  axis  of  the  aircraft.  A  radiation  detector  in  the  focal 
plane  of  the  telescope  converts  the  focused  beam  of  radiation  to  an  electrical  signal.  The 
optical  system’s  field  of  view  first  scans  laterally  across  the  aircraft.  Then,  before  making 
the  next  ground  scan,  it  scans  radiation  references  which  are  internal  to  the  scanner.  By 
the  time  the  next  scan  begins,  the  aircraft  has  moved  forward;  thus  subsequent  line  scans 
build  upon  one  another  to  produce  a  continuous  strip  image  of  the  terrain  beneath  the 
aircraft.  Figure  5.1  shows  the  scanner  looking  directly  downward. 

5.3  -  Preprocessing  Stage 

Because  the  univariate  approximation  chart  of  the  Ozturk  algorithm  is  used  to 
approximate  the  PDF  of  different  regions  the  collected  data  is  therefore  required  to  result 

"  A.  J.  La  Rocca  and  D.  J.  Witte,  “Handbook  of  the  Statistics  of  Various  Terrain  and  Water  (ICE) 
Backgrounds  From  Selected  U.S.  Locations,”  Final  Report.  Environmental  Research  Institute  of  Michigan. 
Contract  No.  N60530-79-R-0036.  January  1980. 
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Figure  5.1  -  Airborne  Scanner  Operation 

in  independent  samples.  Since  there  are  no  tests  available  to  determine  the  statistical 
independence  of  a  set  of  data,  the  best  solution  is  to  ensure  that  the  data  are  spatially 
uncorrelated  in  both  width  and  length.  This  is  done  in  the  preprocessing  stage  of  the 
A’ SCAPE  procedure. 

The  autocovariance  function  has  been  calculated  for  both  sets  of  data  in  the  width 
direction  and  the  length  direction.  Examples  are  shown  in  Figure  5.2  of  the  lag  centered 
at  the  10*  width  pixel  for  all  pixels  in  the  width  direction  for  the  10*  length  pixel  and  in 
Figure  5.3  of  the  lag  centered  at  the  50*  pixel  for  all  pixels  in  the  length  direction  for  the 
10*  width  pixel.  It  has  been  determined  that  in  order  to  obtain  uncorrelated  data, 
decimation  is  necessary  where  only  1  sample  out  of  every  5  samples  is  retained  in  the 
width  direction  and  1  sample  out  of  every  4  samples  is  retained  in  the  length  direction 
reducing  the  scene  size  from  a  high  resolution  383x646  one  to  one  with  a  lower  96x130 
resolution. 
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Note  that  the  decimation  is  necessary  for  the  statistical  procedure  part  of  this  work 
and  is  not  necessary  for  the  mapping  procedure.  Nevertheless,  only  the  data  in  the 
decimated  scene  are  treated  in  this  work  for  the  sake  of  processing  time. 

The  data  processed  in  this  work  and  shown  in  Figure  5.4  is  collected  over  lake 
Michigan  and  consists  of  a  2829  ft  wide  by  1650  ft  long  scene  that  contains  two  major 
regions:  lake  and  land.  The  data  is  8  bits  and  the  resulting  image  size  when  the  pre¬ 
processing  stage  is  completed  is  96x130.  Furthermore,  the  three  dimensional  magnitude 
plot  of  the  scene,  shown  in  Figure  5.5,  reveals  that  the  data  in  the  lake  is  regular  and 
similar  in  magnitudes  whereas  the  data  in  the  land  region  is  irregular  and  bears  a  lot  of 
discretes  especially  near  the  boundary  between  the  land  and  the  lake.  This  latter  fact  is 
due  to  the  non-homogeneity  of  the  land  region  which  contains  a  lot  of  non-homogeneous 
regions  due  to  the  presence  of  trees,  roads,  etc. 

Once  the  pre-processing  stage  is  completed,  the  data  is  forwarded  to  the  mapping 
stage  for  partitioning. 
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Figure  5.4  -  Original  Scene  (Over  Lake  Michigan^ 
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Figure  5.5  -  Three  Dimensional  Magnitude  Plot  of  the  Onpinal  Scene 


5.4  -  Mapping  Procedure  Stage 

After  the  preprocessing  block  has  identified  the  uncorrelated  data  in  the  scene,  the 
mapping  procedure  is  used  to  separate  contiguous  non-homogeneous  regions  with 
different  average  magnitudes. 

When  the  iterative  procedure  starts,  the  diagnosis  procedure  detects  a  discrepancy 
due  to  the  fact  that  the  difference  LPCCP|ncq=8  -  LPCCP|ncq=7  is  less  than  10%  in  all 
iterations.  In  this  case,  and  as  shown  in  Figure  5.6,  LP  (Lake)  and  RPs  (Land)  histograms 
have  a  small  area  of  overlap.  As  summarized  in  Table  4.1,  the  system  searches  for  the 
PLQP  that  results  in  a  minimum  value  LPCCP|f^cQ=g  -  LPCCP|ncq=7- 

After  seven  iterations  of  the  “Identification  of  LP  and  RPs”  stage  the  threshold  is 
accepted  and  the  assessment  passes.  Table  5.1  summarizes  the  values  of  the  different 
parameters  set  in  the  Thresholding/Quantization,  First  and  Second  Correction  stages  as 
described  in  Chapter  4.  Note  that  the  minimum  value  for  LPCCP|ncq=8  -  LPCCP|i.,cq=7  is 
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Figure  5.6  -  Histogram  of  the  Original  Scene 


obtained  for  PLQP=40%  (step  4).  Note  also  that  steps  1  to  5  are  part  of  the  threshold 
approximation  stage  whereas  steps  6  and  7  are  part  of  the  threshold  fine-tuning  stage. 

Figures  5.7  to  5.11  show  the  plots  of  the  QV  (top  right),  CQV  (bottom  left),  and 
CCV  (bottom  right)  scenes  for  steps  1,  2,  3,  4,  and  7,  respectively.  Note  how  the 
quantization  results  more  and  more  in  a  correct  separation  between  the  lake  and  land  as 
PLQP  is  increased.  The  “Identification  of  LP  and  RPs”  stage  is  followed  by  the 
“Detection  of  Patch  Edges”  stage  for  which  the  steps  of  the  last  CCV  (top  left),  SV  (top 
right),  enhanced  edges  volume  (bottom  left),  and  detected  edges  volume  (bottom  right), 
are  shown  in  Figure  5.12  resulting  in  the  segmentation  of  the  scene  into  three  different 
patches  and  their  respective  boundaries.  Note  that  in  addition  to  patches  1  and  2, 
corresponding  respectively  to  land  and  lake,  the  mapping  procedure  detected  a  third 
region,  labeled  3.  This  region  can  be  sighted  in  the  original  scene  of  Figure  5.4  and  is  not 
large.  In  fact  Table  5.2  shows  that  patch  3  contains  only  15  pixels  as  opposed  to  patches  1 
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Step 

LPQP 

NCQ 

NCC 

LPCCP 

LPCCPUcq=8-LPCCP|^cq=7 

(%) 

(%) 

(%) 

1 

10 

8 

1 

24.59 

8.72 

2 

20 

8 

1 

35.82 

6.93 

3 

30 

8 

1 

43.96 

4.91 

B 

40 

8 

1 

46.57 

1.67/ 

5 

50 

8 

1 

56.31 

5.98 

6 

45 

5 

1 

45.14 

- 

7 

45 

5 

2 

45.63 

- 

Table  5.1  -  Setting  of  the  Parameters  in  the  Thresholding/Ouantization,  First  and 


Second  Correction  stages 


Patch  1 

Patch  2 

Patch  3 

Number  of  pixels 

6680 

5337 

15 

Mean 

87.61 

49.47 

34.4 

Variance 

420.20 

4.13 

7.41 

Average  Power  (dB) 

39.08 

33.89 

34.40 

Table  5.2  -  Mapping  Procedure  Assessment 


62 


First  Iteration 


Quantized  Volume  at  10  % 


20  40  60  80  100120 


First  Corrected  Volume  Second  Corrected  Volume 
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Second  Iteration 

LPQP=  20  % 

NCQ  =  8 
NCC  =  1 

LPCCP=  35.82  % 


First  Corrected  Volume 
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Quantized  Volume  at  20  % 
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Second  Corrected  Volume 
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Figure  5.8  -  Step  2  of  the  Identification  of  LP  and  RPs 
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Third  Iteration 

LPQP=  30  % 

NCQ  =  8 
NCC  =  1 

LPCCP=  43.96  % 


First  Corrected  Volume 


20  40  60  80  100120 


Quantized  Volume  at  30  % 
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Second  Corrected  Volume 
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Fourth  Iteration 


Quantized  Volume  at  40  % 
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First  Corrected  Volume  Second  Corrected  Volume 
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Last  Iteration 


Quantized  Volume  at  45  % 
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Figure  5.11  -  Step  7  of  the  Identification  of  LP  and  RPs 
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Last  Corrected  Volume 


Smoothed  Volume 
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Edge-Enhanced  Volume 


Edge-Detected  Volume 
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Figure  5.12  -  Detection  of  Patch  Edges  Stage 
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and  2  which  contain  6680  and  5337  pixels,  respectively.  Furthermore,  the  values  in  Table 
5.2  of  the  variance,  mean  and  average  power  of  patch  3  are  closer  to  those  of  the  lake  than 
those  of  the  land.  This  indieates  that  pateh  3  is  more  related  to  the  lake  than  the  land  and 
may  be  for  example  a  small  body  of  water.  Next,  the  mapping  procedure  is  used  to 
investigate  the  presence  of  any  subpatches  in  every  previously  detected  patch  (i.e., 
patches  labeled  1,  2,  and  3).  In  this  ease,  no  subpatches  are  detected  and  all  three  patehes 
are  declared  as  homogeneous. 

5.5  -  Statistical  Procedure  Stage 

The  statistical  procedure  is  applied  to  every  previously  declared  homogeneous 
patch  in  order  to  separate  further  between  any  existing  contiguous  subpatches  that  may 
have  similar  power  levels  but  different  data  distributions.  Following  the  steps  developed 
in  Chapter  3,  the  data  of  the  scene  is  processed  as  follows: 

Step  1 

Test  pixels  and  their  respective  tiles  (sets  of  100  reference  pixels)  are  selected, 
spread  throughout  the  patch,  to  be  tested  for  Gaussianity  using  the  goodness  of  fit  test  of 
the  Ozturk  algorithm.  Note  that  the  sets  of  100  pixels  are  chosen  to  be  disjoint,  the  closest 
to  and  belonging  to  the  same  patch  as  their  respective  test  pixels.  This  results  in  the  sets 
being  shaped  as  10x10  square  tiles  inside  the  patch  and  tiles  tracking  the  shape  of  the 
boundary  near  the  boundary  of  the  patch.  The  result  of  this  step  is  shovra  in  Figure  5.13, 
where  the  Gaussian  and  non-Gaussian  tiles  are  shaded  in  gray  and  white,  respectively. 
Note  that  many  non-Gaussian  sub-regions  exist  in  patches  1  and  2.  Note  also  that  patch  3 
is  not  processed  by  the  statistical  procedure  due  to  the  fact  that  it  is  composed  of  only  15 
pixels  while  a  minimum  of  100  pixels  are  needed  for  the  Ozturk  algorithm  to  result  in  a 
meaningful  approximation  for  the  distribution  of  the  data. 

Step  2 

Once  the  Gaussian  and  non-Gaussian  regions  are  determined  in  each  patch,  pixels 
with  outliers  are  determined  in  the  non-Gaussian  declared  tiles  (sets  of  100  pixels)  and 
excised. 
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Figure  5.13  -  Scene  After  Step  1  of  the  Statistical  Procedure 


Step  3 

Step  3  of  the  statistical  procedure  proceeds  by  re-checking  for  Gaussianity  the 
previously  non-Gaussian  declared  tiles  when  their  outliers  are  excised.  The  result  is 
shown  in  Figure  5.14.  Note  that  even  though  a  lot  of  non-Gaussian  previously  declared 
tiles  are  now  declared  Gaussian,  non-Gaussian  tiles  still  exist  particularly  in  the  land 
region  near  the  boundary  with  the  lake.  This  is  because,  as  stated  previously,  many  non- 
homogeneous  regions  due  mainly  to  trees  and  roads  exist  in  that  region  causing  discretes 
to  appear  in  the  data  thus  making  it  non-Gaussian.  Figure  5.15  shows  the  location  of 
outliers  (black  dots).  Comparing  Figures  5.15  and  5.4,  note  that  outliers  tend  to  have  a 
physical  significance  and  might  represent  string-like  patches  such  as  road-ways. 

The  results  of  the  statistical  procedure  (steps  1  to  3)  are  summarized  in  Table  5.3. 
Note  that: 


70 


Figure  5.14  -  Scene  After  Step  3  of  the  Statistical  Procedure 


20  40  60  80  100  120 


Figure  5.15  -  Scene  After  Step  3  of  the  Statistical  Procedure  and 

Location  of  Outliers 
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(1)  the  outliers  occupy  a  small  percentage  of  the  scene  (1.63%), 

(2)  Gaussian  regions  occupy  the  largest  part  of  the  scene,  and 

(3)  non-Gaussian  regions  exist  and  occupy  about  16%  of  the  scene. 


Percentage  of 
Gaussian 
Pixels  (%) 

Percentage  of 
Non-Gaussian 
Pixels  (%) 

Percentage 
of  Outliers 
(%) 

Tiles  with  present 
outliers 

61.81 

38.19 

- 

Tiles  with  excised 
outliers 

82.43 

15.94 

1.63 

Table  5.3  -  Percentages  of  Gaussian  and  Non-Gaussian  regions 


Step  4 

In  this  step,  the  following  procedure  is  followed, 

1  -  non-Gaussian  tiles  are  numbered  as  shown  in  Figure  5.16. 

2  -  contiguous  numbered  tiles  are  then  grouped  in  sets  as  shown  in  column  1  of 
Table  5.4. 

3  -  Using  the  Ozturk  algorithm,  the  (U,V)  coordinates  of  the  locus  end  point  is 
obtained  for  every  tile. 

4  -  The  average  (Uav,Vav)  coordinates  is  computed  for  every  set  of  contiguous 
non-Gaussian  tiles  and  the  average  approximating  PDF  corresponding  to  the  pair 
(Uav,Vav)  is  obtained.  The  PDF  type  along  with  the  first  and  second  shape  parameters 
are  shown  in  columns  3,  4,  and  5  of  Table  5.4,  respectively  for  every  set  of  tiles.  Note 
that  a  set  can  consist  of  a  single  tile  if  the  tile  is  not  contiguous  to  any  other  non-Gaussian 
tiles. 
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Tile 

Fitting  under 

PDF  Type 

1st  Shape 

2nd  Shape 

Number 

same  ellipse 

Parameter 

Parameter 

1 

- 

Beta 

6.0 

1.6 

2 

- 

Beta 

1.0 

3.2 

3 

- 

SU-J 

1.2 

-0.2 

IB 

Yes 

Beta 

1.0 

3.2 

6,  7,8 

No 

Beta 

2.0 

3.2 

9  to  16 

Not  all 

Beta 

1.0 

3.2 

17 

- 

Beta 

1.0 

0.8 

18 

- 

Beta 

2.0 

0.8 

19 

- 

Beta 

2.0 

1.6 

20,21 

Yes 

Beta 

2.0 

1.6 

22 

- 

SU-J 

1.2 

-0.2 

23 

- 

Beta 

1.0 

1.6 

Table  5.4  -  PDFs  to  approximate  the  Different  Subpatches 


5  -  A  check  is  made  to  ensure  that  all  (U,V)  pairs  of  every  set  fit  withing  the  same 
confidence  ellipse  as  (Uav,Vav).  The  result  of  this  check  is  reported  in  column  2  of  Table 
5.4.  Note  that  two  sets,  consisting  of  tiles  numbered  6  to  8  and  9  to  16,  do  not  pass  the 
check. 

6  -  For  those  sets  that  do  not  pass  the  check  above,  the  tiles  are  regrouped  so  that 
all  (U,V)  coordinates  for  each  group  of  tiles  can  fit  within  the  same  ellipse.  Each  group  of 
contiguous  tiles  form  then  a  subpatch.  Table  5.5  shows  the  result  of  this  step.  Note  that 
column  1  displays  the  sets  of  contiguous  tiles  whereas  column  2  allocates  numbers  to  the 
subpatches  obtained  by  the  statistical  procedure.  Note  also  that  each  subpatch  has  its  PDF 
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Tile 

Sub-patch 

Fitting  under 

PDF  Type 

1st  Shape 

2nd  Shape 

Number 

Number 

same  ellipse 

Parameter 

Parameter 

1 

1 

- 

Beta 

6.0 

1.6 

2 

2 

- 

Beta 

1.0 

3.2 

3 

3 

- 

SU-J 

1.2 

-0.2 

4,5 

4 

Yes 

Beta 

1.0 

3.2 

6 

5 

- 

Beta 

1.0 

0.8 

7 

6 

- 

Beta 

1.0 

1.6 

8 

7 

- 

Lognormal 

0.6 

N/A 

9  to  15 

8 

Yes 

Beta 

1.0 

3.2 

16 

9 

Beta 

0.6 

0.8 

17 

10 

- 

Beta 

1.0 

0.8 

18 

11 

- 

Beta 

2.0 

0.8 

19 

12 

- 

Beta 

2.0 

1.6 

20,21 

13 

Yes 

Beta 

2.0 

1.6 

22 

14 

- 

SU-J 

1.2 

-0.2 

23 

15 

- 

Beta 

1.0 

1.6 

Table  5.5  -  PDFs  to  approximate  the  Different  Subpatches 


approximated  as  shown  in  columns  4  to  6  of  Table  5.5. 

At  the  end  of  the  statistical  procedure,  subpatches  within  the  lake  and  land  regions 
are  isolated.  Outliers  which  represent  tiny  patches  such  as  road  ways  are  located  and  for 
every  pixel  in  the  scene,  the  PDF  to  approximate  the  data  is  readily  known. 
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5.6  -  Indexing  Stage 

The  statistical  procedure  stage  is  followed  by  the  indexing  stage  in  which  every 
pixel  in  the  scene  is  assigned  a  set  of  descriptors  to  indicate,  for  each  pixel,  (1)  its 
identification  number  which  is  the  same  as  the  identification  number  of  the  homogeneous 
patch  or  subpatch  to  which  the  pixel  belongs,  (2)  whether  the  pixel  is  an  inner  pixel  or 
edge  pixel,  (3)  its  type  (Gaussian,  non-Gaussian,  or  outlier),  (4)  for  a  non-Gaussian  pixel, 
its  best  approximating  PDF. 

Up  to  this  point,  patches  1,  2,  and  3  have  been  identified,  outliers  have  been 
located  in  patches  1  and  2.  Also,  subpatches  have  been  defined  within  every  one  of 
patches  1  and  2  and  PDFs  to  approximate  the  different  non-Gaussian  subpatches  are 
determined  using  the  application  chart  mode  of  the  Ozturk  algorithm.  All  of  this 
information  is  readily  available  and  stored  in  the  indexing  stage. 

5.7  -  Conclusion 

The  processing  of  the  real  IR  data  by  A’ SCAPE  has  resulted  during  the  mapping 
procedure  in  the  partitioning  of  the  scene  into  three  main  patches,  lake,  land,  and  a 
subpatch  within  the  land.  Then,  the  statistical  procedure  partitioned  further  the  lake  and 
land  into  main  Gaussian  patches  and  a  total  of  fifteen  non-Gaussian  subpatches.  Also, 
outliers  representing  tiny  patches  such  as  roads  have  been  located  and  approximating 
PDFs  have  been  determined  for  the  non-Gaussian  subpatches.  Finally,  if  A’ SCAPE  is  to 
be  followed  by  a  detection  stage,  all  information  pertaining  to  what  type  of  detector 
should  be  used  and  what  the  parameters  should  be  is  readily  available. 
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6 


A’SCAPE  Demo  Package 


6.1  -  Introduction 

A  demo  package  is  built  in  Matlab  which  describes  in  detail  the  different  stages  of 
A’SCAPE.  The  package  has  a  friendly  mouse  driven  graphical  user  interface  (GUI)  and 
consists  of  two  main  sections.  The  first  section  presents  the  detailed  steps  of  the  real  IR 
data  example  of  Chapter  5.  The  second  section  is  subdivided  into  two  subsections  where 
in  the  first  one  a  set  of  examples  is  presented  which  illustrate  the  need  Tor  A’SCAPE, 
whereas  in  the  second,  a  detailed  description  is  given  for  every  stage  of  A’SCAPE.  The 
views  can  be  displayed  manually  or  in  an  automated  way  as  a  slide-show.  A  movie  is 
included  which  shows  the  steps  in  the  mapping  procedure.  The  next  sub-sections  describe 
in  detail  the  different  parts  of  the  demo  package. 

As  shown  in  Figure  6.1,  note  that  the  A’SCAPE  demo  consists  of  three  windows 
referred  to  as  the  menu  window  (top  right),  plot  window(top  left),  and  command  window 
(bottom  left).  These  are  described  next. 

6.2  -  Menu  Window 

The  menu  window  consists  of  a  main  menu  and  sub-menus  formed  each  of  a  set 
of  vertically  aligned  push  buttons.  Each  button  when  pressed  either  generates  another 
menu  or  plots  in  the  Plot  Window.  The  first  menu  to  appear  in  the  demo  is  referred  to  as 
the  Main  Menu.  As  shown  in  Figure  6.2,  the  first  three  buttons  of  the  main  menu  are 
labeled  Original  Scene,  Mapping  Procedure,  and  Statistical  Procedure,  respectively. 
These  buttons  run  the  example  of  the  real  IR  data  scene  treated  in  chapter  5.  The  fourth 
button,  labeled  A’SCAPE  Tutorial,  offers  a  tutorial  on  A’SCAPE.  The  fifth  button, 
labeled  About  A’SCAPE,  gives  credit  to  the  author  and  sponsors  of  the  work.  Finally,  the 
button  labeled  end  stops  the  demo  and  closes  all  windows. 

All  of  the  first  four  buttons  generate  sub-menus  labeled  as  summarized  in  Table 

6.1. 
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Menu  Window  (top  right).  Plot  Window  (top  left).  Command  Window  (bottom  left 


Figure  6.2  -  A’SCAPE  Demo:  Main  Menu 
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Main  Menu 
Button  Number 

Effects 

1 

Sub-menu  of  Figure  6.3  labeled  Original  Scene  Menu 

2 

Sub-menu  of  Figure  6.4  labeled  Mapping  Procedure  Menu 

3 

Sub-menu  of  Figure  6.5  labeled  Statistical  Procedure  Menu 

4 

Sub-menu  of  Figure  6.6  labeled  A’ SCAPE  Tutorial  Menu 

Table  6.1  -  Effects  of  Buttons  1  to  4  of  the  Main  Menu 

Effects  of  the  buttons  in  each  of  the  sub-menus,  generated  by  buttons  1  to  4  of  the 
main  menu  and  shown  in  Figures  6.3  through  6.5,  are  summarized  in  tables  6.2  to  6.5, 
respectively. 


Figure  6.3  -  Original  Scene  Menu 
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Figure  6.4  -  Maf 


Procedure  Menu 


Figure  6.5  -  Statistical  Procedure  Menu 
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Original  Scene 
Menu 

Button  Label 

Effects 

Image 

Displays  the  image  of  the  original  scene 

Contour 

Displays  the  contour  plot  of  the  original  scene 

Magnitude 

Displays  the  3  dimensional  plot  of  the  original  scene 

Histogram 

Displays  the  histogram  of  the  original  scene 

All  of  the  Above 

Displays  all  of  the  above  in  a  single  window 

Main  Menu 

Closes  Original  Scene  Menu  and  opens  Main  Menu 

Table  6.2  -  Effects  of  Buttons  in  the  Original  Scene  Menu 
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Mapping  Procedure 
Menu 

Button  Label 

Effects 

Quantization 

Displays  the  quantization  steps 

Correction 

Displays  the  correction  steps 

Edge  Enhancement 

Displays  the  edge  enhancement  steps 

Last  Iteration 

Displays  the  steps  in  the  last  iteration 

Mapping  Assessment 

Displays  the  steps  in  the  assessment  of  the  Mapping 
Procedure 

Sub-Patch 

Investigation 

Displays  the  steps  of  the  investigation  of  sub-patches 

All  of  the  Above 

Displays  all  of  the  above  in  a  sequential  manner 

Summary 

Displays  all  of  the  above  in  a  compact  manner 

Movie  of  Mapping 

Displays  the  mapping  procedure  as  a  movie 

Main  Menu 

Closes  Mapping  Procedure  Menu  and  opens  Main 
Menu 

Table  6.3  -  Effects  of  Buttons  in  the  Mapping  Procedure  Menu 


Statistical  Procedure 
Meuu 

Button  Label 

Effects 

Non-Gaussian 

Regions 

Displays  the  search  for  Gaussian  and  Non-Gaussian 
region  steps 

Region  PDF 

Approximation 

Displays  the  search  for  PDFs  approximating  Non- 
Gaussian  region 

Assessment 

Assessment  of  the  Statistical  Procedure 

All  of  the  Above 

Displays  all  of  the  above  in  a  sequential  manner 

Summary 

Displays  all  of  the  above  in  a  compact  marmer 

Main  Menu 

Closes  Statistical  Procedure  Menu  and  opens  Main 
Menu 

Table  6.4  -  Effects  of  Buttons  in  the  Mapping  Procedure  Menu 
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A’SCAPE  Tutorial 
Menu 

Button  Label 

Effects 

Examples 

Displays  examples  to  justify  the  need  for  A’SCAPE 

A’SCAPE 

Displays  block  diagram  of  A’SCAPE 

Mapping  Procedure 

Displays  block  diagram  of  the  Mapping  Procedure 

Statistical  Procedure 

Displays  block  diagram  of  the  Statistical  Procedure 

Approximation  Chart 

Displays  an  example  of  an  approximation  chart 

IPUS 

Displays  a  summary  of  the  expert  system  main 
definitions 

All  of  the  Above 

Displays  all  of  the  above  in  a  sequential  manner 

Main  Menu 

Closes  A’SCAPE  Tutorial  Menu  and  opens  Main 
Menu 

Table  6.5  -  Effects  of  Buttons  in  the  Mapping  Procedure  Menu 


6.3  -  Plot  Window 

All  plots  displayed  by  the  different  buttons  are  presented  in  the  Plot  Window.  The 
plots  are  either  full  sized  one  plot  at  a  time  or  compact  four  plots  at  a  time. 

6.4  -  Command  Window 

When  there  is  more  than  one  plot  to  be  displayed  by  a  single  button  in  any  of  the 
menus,  buttons  appear  in  the  Command  Window  as  shown  in  Figure  6.7  which  consist  of 
four  push  buttons  and  two  sliders.  Effects  of  the  buttons  are  described  in  Table  6.6. 


Figure  6.7  -  Buttons  in  the  Command  Window 
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Command  Window 
Button  Label 

Effects 

Step 

Displays  frame  by  frame  at  the  request  of  the  user 

Auto 

Displays  the  frames  automatically  with  a  pause 
between  the  frames  defined  by  the  Pause  slider 

Stop 

Stops  the  automatic  display  of  the  different  frames 

Done 

Clears  the  command  window  and  returns  command  to 
the  menu  window 

Pause  (slider) 

Shows  the  pause  in  seconds  between  the  frames  being 
displayed  automatically.  Also  enables  choice  of  the 
pause  in  seconds  when  the  slider  is  moved 

Frame  (slider) 

Shows  the  number  of  the  frame  being  displayed.  Also, 
enables  display  of  any  desired  frame  when  the  slider 
is  moved 

Table  6.6  -  Effects  of  Buttons  in  the  Mapping  Procedure  Menu 
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Part  II 

Performance  Analysis  of  Multichannel  Parametric 
Detection  Algorithms 
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Introduction 

In  many  signal  processing  techniques,  the  number  of  signals,  say  q,  is  assumed  to  be 
known  a  priori.  This  is  needed  for  further  processing  of  the  received  data  vector.  However, 
in  practice,  q  is  unknown  and  therefore  needs  to  be  estimated.  This  problem  is  often  referred 
to  in  the  literature  as  a  "Detection  Problem".  For  simplicity,  we  review  the  basis  of  the 
subspace  approaches.  This  is  needed  for  a  good  understanding  of  our  work.  Let  the  number 
of  elements  in  an  array  of  sensors  be  denoted  by  m,  where  m  is  assumed  to  be  greater  than  q. 
The  received  signal  vector  can  be  expressed  as 

X=AS+N  (1-1) 

where  X  and  N  are  (mxl)  complex  vectors,  S  is  a  (qxl)  complex  vector  and  A  is  an  (mxq) 
complex  matrix  containing  information  about  the  signals  of  interest.  S  refers  to  the  random 
signal  vector  while  N  is  a  noise  vector.  All  signal  and  noise  wave  forms  are  assumed  to  be 
zero-mean  random  narrow  band  processes.  Moreover,  for  simplicity,  the  additive  noise  vector 
is  assumed  to  have  a  diagonal  covariance  matrix  given  by  where  I^  is  the  (mxm) 

identity  matrix.  The  case  of  correlated  noise  is  considered  m  section  2.2.2.  We  also  assume 
that  the  signal  and  noise  vector  are  statistically  independent.  Note  that  the  analysis  considers 
the  spatial  domain  only. 

It  follows  from  the  signal  model  represented  by  Equation  (1-1)  and  the  above 
assumptions  that  the  covariance  matrix  of  X  is  given  by 

Rx  =E[X«  X'']  =  ASA”  =Rs  +R^  =  Rs  (1-2) 

where  Rg  is  the  covariance  matrix  of  the  signal  vector  S  and  H  denotes  the  complex  conjugate 
transpose  operation.  In  an  eigen-decomposition,  note  that  the  rank  of  the  (mxm)  matrix  R^  is 
q.  Therefore  the  smallest  eigenvalues  of  R^  are  identical  and  are  equal  to  .  If  the 
eigenvalues  of  R^  are  ordered  in  descending  order  with  respect  to  their  magnitudes,  such  that 
X.]  >^2  >"->X^  ,  it  follows  that  the  q  largest  eigenvalues  of  R^  have  magnitudes 

greater  than  a^.  Specifically, 
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>cl  ;i=l,2, q 


(1-3) 


and 

A,;  =  ;  i  =  (q+1),  (q+2), m.  (1-4) 

If  the  exact  eigenvalues  of  the  covariance  matrix  of  the  received  signal  were 
available  to  us,  the  solution  to  the  problem  of  detection  (estimating  q)  would  be 
straightforward  in  that  we  would  simply  have  to  count  the  multiplicity  of  the  smallest 
eigenvalue.  In  practice,  however,  the  covariance  matrix  and  its  eigenvalues  are  unknown 
and  they  have  to  be  estimated  from  a  finite  set  of  data  samples.  Due  to  the  noisy  nature  of 
these  estimates,  we  find  that  the  exact  multiphcity  rule  does  not  hold.  Hence,  this  simple  rule 
for  determining  the  number  of  targets  has  to  be  dropped  in  favor  of  a  more  robust  approach. 

The  first  technique  to  be  used  in  the  detection  problem  was  based  on  a  statistical 
analysis  of  the  received  data.  Various  criteria  are  described  in  the  hterature  for  determining 
the  number  of  targets.  The  two  most  commonly  used  are  the  Akaike  information  criterion 

(AIC)  introduced  by  Akaike  1,  which  was  first  introduced  in  Direction  Finding  (DF)  problems 
by  Wax  et  al^,  and  the  Minimum  Descriptive  Length  (MDL)  criterion  of  Shwartz^  and 
Rissanen^. 

In  both  approaches,  the  signals  and  m  additive  noise  wave  forms  are  assumed  to  be 
ergodic  Gaussian  processes.  Let.  R^  denote  the  Maximum  Likelihood  (ML)  estimate  of  the 

covariance  matrix  R^.  Also,  let  ;k  =  l,2,-",m,  1,  >4  be  the  eigenvalues  of 

Rj^ .  The  cost  function  for  both  the  AIC  and  MDL  criterion  takes  on  the  form 


^  H.  Akaike,  "A  New  Nook  at  the  Statistical  Model  Identification,"  IEEE  trans.  Automat.  Contr.,  Vol.  AC  - 
14,  pp.  716-723,  Dec.  1974. 

2  M.  Wax  and  T.  Kailath,  "Detection  of  Signals  by  Information  Theoretic  Criteria,"  IEEE  Trans.  Acoust., 
Speech.  Signal  Processing.  Vol.  ASSP-33,  pp.  387-392,  Apr.  1985. 

2  G.  Schwartz,  "Estimating  the  Dimension  of  a  Model."  Anna.  Stat..  Vol.  6,  No.  22,  pp.  461-464,  1978. 

J.  Rissanen,  "A  Universal  Prior  for  Integers  and  Estimation  by  Minimum  Descriptive  Length,"  Ann.  Stat., 
Vol.  11,  No.  2,  pp.  466-471,  1983. 
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J(k)  =  L(k)+p(k) 


(1-5) 


where  the  log-likelihood  function,  L(k)  is  given  by 


L(k)=N-Ln 


1 


I  m-k 


(m-k), 5''; 


ni, 

i=K=l 


(1-6) 


with  N  being  the  number  of  time  samples  (snapshots)  used  in  evaluating  ,  and  the  penalty 
function,  p(k),  is  given  by 

p(k)=a(N)[k(2m-k)],  (1-7) 

with  a(N),  the  penalty  coefficient,  being  a  function  of  the  number  of  samples. 

The  choice  of  a(N)  depends  on  the  criterion  used.  For  the  AIC,  a(N)  =  l  whereas  for 
the  MDL  criterion,  a(N)  =  (l/2)ln(N) .  The  number  of  signals  q  is  then  determined  as  that 
value  of  k  for  which  J(k)  is  minimized. 

Let  Hq  denote  the  hypothesis  that  the  true  number  of  targets  is  q.  Then,  the 
probability  of  underestimating  and  overestimating  the  number  of  targets,  given  ,  are 


PM=p(q<a/Hj 

(1-8) 

Pp=p(q>q/H,) 

(1-9) 

P^andPp  are  called  the  miss  and  false  alarm  probabilities,  respectively.  It  has  been 

shown^  that  the  probability  of  missing  for  both  the  AIC  and  MDL  criterion  can  be  reduced  to 
zero  by  either  increasing  the  signal-to-noise  ratio  (SNR)  or  by  increasing  the  number  of 
samples  N.  However,  for  both  these  parameters,  the  decrease  is  much  faster  for  AIC  than  for 
MDL.  This  indicates  that  the  AIC  is  more  efficient  in  reducing  the  probabihty  of  missing  than 


5  Q.  Zhang,  K.  Wang,  Y.  Yin,  J.  Reilly,  "Statistical  Analysis  of  the  Performance  of  Information  Theoretic 
Criteria  in  the  Detection  of  the  Number  of  Signals  in  Array  Processing,"  IEEE  Trans.  Acoust.,  Speech,  Signal 
Processing.  Vol.  ASSP-37,  pp.  1557-1567,  Oct.  1989. 
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is  the  MDL  criterion.  On  the  other  hand,  for  a  moderate  number  of  samples  the  probabihty  of 
false  alarm  using  the  MDL  criterion  approaches  zero  with  increasing  values  of  N,  whereas  that 
for  the  AIC  remains  constant.  This  shows  that  the  MDL  criterion  is  more  efficient  in  reducing 
the  probability  of  false  alarm  than  is  the  AIC.  Under  low  SNR’s,  both  the  AIC  and  the  MDL 
necessitate  a  large  number  of  samples,  which  may  be  very  hard  to  meet  in  practice. 

Research  is  needed  to  determine  whether  the  performance  of  the  information  theoretic 
criteria  can  be  improved  for  specific  applications  by  selecting  other  functions  for  a(N)  and/or 
log-likelihood  functions  L(k). 

Another  technique^,  based  on  the  AIC  and  the  MDL,  has  been  applied  on  the  singular 
values  of  the  covariance  matrix  instead  of  the  eigenvalues  of  the  covariance  matrix.  The 
reason  is  that  there  is  a  direct  relationship  between  the  eigenvalues  of  the  covariance  matrix 
and  its  singular  values.  The  advantage  of  this  technique  over  the  previous  one  is  that  it 
benefits  from  the  numerical  efficiency  and  stability  of  the  SVD  technique. 

Recently,  a  new  procedure^  for  the  detection  of  the  number  of  targets  was  introduced, 
within  the  framework  of  model  selection,  and  established  the  strong  consistency  of  the 
procedure.  This  procedure  is  referred  to  as  the  Efficient  Detection  Criterion  (EDC) 
procedure.  According  to  this  procedure,  q  is  estimated  with  q  where 

EDC(q)  =  min{EDC(0),  .  .  .  ,  EDC(m-l)}  (1-10) 

where 

EDC(k)  =  -2Ln(k)  +  v(k,m)C(N)  (1-11) 

and  C(N)  satisfies  the  following  conditions: 

y  lim[c(N)/N]=0 
N  — ^  °° 

6  V.  Shamirian  and  S..  B.  Kesler,  “Detection  of  Signals  by  SVD-Based  Information  Theoretic  Criteria,”  Proa 
30-th  Midwest  Symposium  on  Circuits  and  Systems,  pp.  567-570,  Syracuse,  NY  1987. 

L.  C.  Zhao  et  al.,  “Remarks  on  Certain  Criteria  for  Detection  of  Number  of  Signals,”  IEEE  Trans.  ASSP, 
Vol.  ASSP-35,  No.  11,  pp.  129-132,  Eeb.  1987. 
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(1-13) 


lim{c(N)/Ln[Ln(N)]}=0 

N  — >  oo 

A  very  general  and  elegant  method^  to  estimate  the  number  of  targets  was  recently 
introduced.  This  method  depends  on  choosing  a  function  f  defined  on  aU  finite  sequences, 
called  r-regular  functions.  The  previous  methods  used  only  1 -regular  functions  (r  =  1).  It  was 
shown  that  actually  this  is  the  worst  case.  It  was  eventually  shown  that  as  r  increases,  the 
convergence  is  even  faster. 

From  the  above  discussion,  we  notice  that  most  techniques  are  non-parametric.  It  was 
our  belief  that  parametric  techniques  are  better  suited  for  this  type  of  detection  since  they  have 
more  degrees  of  freedom  and  eventually  will  prove  very  efficient. 

Our  interest  in  the  parametric  techniques  and  LeCadre's  technique^  in  particular  is  two 
fold.  First,  we  evaluated  and  studied  its  performance  against  several  parameters  such  as  order 
determination  of  the  AR  model,  number  of  incoming  signals  (targets),  signal-to-noise  ratios, 
and  the  presence  or  absence  of  clutter.  We  also  studied  the  effects  of  a  limited  data  sequence, 
computational  complexity,  algorithm  convergence  and  implementation  feasibility.  This  step 
could  be  referred  to  as  a  complete  assessment  of  the  proposed  technique.  In  the  second 
phase,  we  mainly  restrict  ourselves  to  the  software  development  of  appropriate  detection  and 
estimation  algorithms  for  insertion  in  the  Rome  Lab  Space  Time  Signal  Processing  (RLSTAP) 
software. 


^  Y.  Q.  Yin  and  P.  K.  Krishnaiah,  “On  Some  Non-Parametric  Methods  for  Detection  of  the  Number  of 
Signals,”  IEEE  Trans.  ASSP.  Vol.  ASSP-35,  No.  11,  pp.  1523-1538,  Nov.  1988. 

^  J.  P.  LeCadre,  “Parametric  Methods  for  Spatial  Signal  Processing  in  the  Presence  of  Unknown  Colored 
Noise  Fields,”  IEEE  Trans.  ASPP.  Vol.  37,  No.  7,  pp.  965-983,  July  1989. 
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2 


LeCadre’s  Technique 

In  this  section,  we  will  derive  the  detection  algorithm  defined  by  LeCadre  and  provide 
insights  into  the  models  used  in  the  derivation. 

2.1  -  Introduction 

Let  Xj  be  the  i-th  stationary,  m-dimensional  random  vector  constituted  by  the  narrow- 
band  outputs  of  an  array  sampling  a  homogeneous  random  field,  where  m  is  the  number  of 
sensors.  Xi  can  be  written  as 

Xi=[xi(i)  x^Ci)  ...  x,„(i)f.  (2-1) 

Let  be  the  covariance  matrix  of  the  outputs  vector  Xj .  can  be  written  as 

Rx=E[Xi*Xr]-  (2-2) 

It  can  easily  be  shown  that  R^  can  be  expressed  as 

Rx=Rs+Rn’  (2-3) 

where  Rj  and  R^^  are  the  covariance  matrices  of  the  sources  and  noise,  respectively. 

The  objective  is  to  obtain  an  accurate  estimation  of  R^,  based  on  the  availability  of  an 
estimate  of  the  covariance  matrix  Rj, .  The  method  relies  on  3  facts: 

1.  The  likelihood  functional  may  be  expressed  as  a  function  of  the  eigenvalues  of  the 
whitened  covariance  matrix  of  the  outputs. 

2.  The  inverse  of  the  noise  covariance  matrix  admits  an  explicit  formulation  in 
terms  of  the  AR(MA)  coefficients  of  the  noise  model. 

3.  The  derivatives  of  the  likelihood  functional  may  be  computed  using  (1)  and  (2)  and 
classical  results  for  perturbations  of  eigenvalues. 
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2.2  -  Calculation  of  the  Likelihood  Functional 

The  derivation  of  the  likelihood  functional  depends  upon  which  type  of  interference 
(i.e.;  thermal  white  noise  only  or  correlated  clutter  plus  additive  white  noise)  is  being  used. 
The  first  case  will  be  devoted  to  the  white  noise  interference  and  the  second  case  will  be 
devoted  to  the  non-white  case. 


2.2.1  -  White  Noise  case 

In  this  case,  we  assume  that  is  of  the  form  where  is  a  positive 

constant  and  is  the  (mxm)  identity  matrix.  Moreover,  assume  that  the  number  of  sources 
is  denoted  by  q.  Let  {X,  ,X2  be  a  sequence  of  N  independent  complex  Gaussian 

vectors  with  covariance  matrix  R.  Note  that  R  is  given  by 


R.=;!rix,‘xr. 


where  H  denotes  complex  conjugate  transpose. 

Given  Rg  and  ,  the  pdf  of  any  spatial  vector  snapshot  Xj  is  given  by 

f.ixAH.Rs)=^^exp{-xrR;'x,},  (2-5) 

where  independence  was  assumed  from  sensor  to  sensor  and  over  the  sequence  of  N  samples. 
The  conditional  density  of  the  received  sequence  is  then  given  by 


>^s)“rT^Xi  (ii/^N  ’^s) 


W""  [de.(R.)]' 


expj-X^fRx’Xi 


To  simplify  this  expression,  we  study  the  following  example.  Consider  the  case  where 


X2]^  and  let  R,^'  = 


N 

Note  that  is  equal  to 

i=l 
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i=l 


i=l 


r  * 

*  1 

[a 

b‘ 

- 1 

>< 

_ 1 

Fii 

XizJ 

I - 

o 

d 

1 - 

to 

1 _ 

+bXi2x;,  +cXiiX*2  +dXi2X*2 .  (2-7) 


i=l 


On  the  other  hand,  note  that 

N  N 

NR,=ix,.xr=x 

i=l  i=l 

Therefore  the  product  NR  ^  R  is  given  by 


- 1 

' 

_ I 

[x-l  Xi2]  =  X 

Xi,Xi,  x„x;2 
♦  * 

L^i2j 

i=l 

LXi2Xi,  Xi2Xi2j 

nr.r;'=S 


i=l 


Xn^ii  Xi,Xi2 

Xi2Xn  Xi2X;2 


a  b 
c  d 


=I 


i=l 


aXi,Xi, -l-cXjiXi2  bXi,x*, -l-dXi,Xi2 

aXj2X;,  +  CXi2x‘2  bx.2X;,  +  dXi2X;2  J 


(2-8) 


(2-9) 


Recall  that  the  trace  of  a  matrix  A  is  defined  as  the  sum  of  its  diagonal  elements;  i.e., 
tr(A)  =  ^aii  .  Note  in  this  case  that  the  trace  of  NR^^R"'  is  given  by 


iN 

tr(NR,R;’)=Ntr(R,R;')=XaXiiXi, -Hcx„x*2+bXi2X*, +dXi2Xj2  ■  (2-10) 


Note  that  Equations  (2-7)  and  (2-10)  are  identical.  Hence,  Equation  (2.2. 1-3)  can  be  written 
as 


fx(2iAN>Rs)=7^i— 7^7;^exp{-Ntr(R,R;’)}.  (2-11) 

in)  [det(Rj] 

The  next  step  is  to  find  the  value  of  and  Rg  which  maximize  the  likelihood 
functional.  Note  that  Ln[f ^  (x  /  iV  ^ ,  Rs )]  is  given  by 

Ln[fx  (x/^^N  ’  Rs  )]= “  N  jmLn(7t)  -I-  Ln[det(R  J]  -f-  tr(R,,R;' )} .  (2-12) 

It  can  be  shown  that  maximizing  this  function  with  respect  to  results  in  the  following 
equation 


‘0  T.  W.  Anderson,  “Asymptotic  Theory  for  Principal  Component  Analysis,”  Ann.  J.  Math.  Stat.,  Vol.  34,  pp. 
122-148,  1963. 
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m 

maxjLnj^fx  (x/^n  > )  j  =  Nj  -  mLn(7i)  -  Ln[det(R^ )]  +  ^ Ln(Xi )  -(m-q)Ln 


i=q+l 


Let  C  be  the  following  constant 


C  =  -  mLn(7t)  -  Ln[det(R  ^ )] , 


-I  m 

—  X^I 

m-qi=q+i 

(2-13) 

(2-14) 


then,  the  above  equation  can  be  expressed  as 


(1  m  ,  ^  » 

C  +  (m-q)— — Ln  n(^i) 

Qj  _i=q+l'' 


-  (m  -  q)Ln 


1 


(2-15) 


This  last  equation  can  also  be  written  as 


max  {^nj^fx  (x/A<[q  .Rj  jj|— N 


C  -h  (m  -  q)Ln 


n(M 

i=q+l 


(/m-q 


-  (m  -  q)Ln 


1 


m-qi^, 


(2-16) 

Therefore,  it  can  easily  be  seen  that  maximizing  the  likelihood  amounts  to  maximizing  the 
functional 


Lq  (^N  >  Rs  )= (m  -  q)Ln[ge(q)]  -  (m  -  q)Ln[ar(q)] ,  (2-17) 

where  ge(q)  and  ar(q)  are  given  by 


and 


(2-18) 


1  ^ 

ar(q)= - • 

m-qi=q+, 


(2-19) 
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ge(q)  and  ar(q)  are  referred  to  as  the  geometric  and  arithmetic  means,  respectively.  Note  also 
that  in  the  previous  derivation,  we  have  assumed  that  the  eigenvalues  of  are  ordered;  i.e., 

(2-20) 

We  now  consider  the  case  of  unknown  correlated  noise. 


2.2.2  -  Unknown  Correlated  Noise 

In  this  case,  the  likelihood  functional  is  given  by 

L,(R„,R,)=-Ln[det(Rj]-tr(R;'R.).  (2-21) 

Note  that  is  a  positive  definite  matrix,  therefore  it  admits  a  Cholesky  decomposition  in  the 
form  of  two  triangular  factors  (matrices).  Let  R^  be  given  by  R,^  =L.L”  .  Consider  then  the 
transformed  matrix  of  the  exact  covariances  (after  application  of  L);  i.e., 

R^  =  L-'RL-"  =  L-'  (Rs  +  XnRn  )L'”  =  L"’  (Rg  +  )l-"  ,  (2-22) 

which  can  also  be  written  as 

R^  =  L-'RjL-”  +  A,nL-'LL"L-”  =  L-’RjL-”  (2-23) 

Furthermore,  let  be  the  transformed  source’s  matrix;  i.e.,  =L“’RsL““.  Then  the 

factor  tr(R“'R,^)  becomes 

tr{R;'Rj  =  tr[(L-"R;'L-')(LR.L")]=tr(L-"R;'R.L“).  (2-24) 

For  simpMcity,  let  A  and  B  be  the  following  matrices,  A  =  L“”R”'  and  B  =  R„L“.  Also, 
recall  that  given  two  matrices  A  and  B,  it  is  known  that  tr(AB)=  tr(BA).  Therefore  Equation 
(2-24)  becomes 

tr(R;'R.)  =  tr(L-“R;'R.L”)=tr(R,L"L-"R;')=tr(R.R;')=tr(R;‘R,).  (2-25) 

Note  also  that  the  determinant  of  Rx,  det(Rx),  can  be  written  as 
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det(R  J  =  det(LR^L" )  =  det(L)  •  det(L” )  •  det(R„ )  =  det(R^ )  •  det(R„ ) ,  (2-26) 

since  R^  =L.L^  and  det(R^ )  =  det(L)  •  det(L“ ) .  Therefore  Ln[det(Rx)]  is  given  by 

Ln[det(R^ )]  =  Ln[det(R^ )]  -I-  Ln[det(R^ )] .  (2-27) 

Thus,  Equation  (2-21)  becomes 

Lq  (Rn  .  Rs  )=-  Ln[det(R„  )]  -  tr(R;'  R„ )  -  Ln[det(RN )].  (2-28) 

Let  (q)  be  the  quantity 

Lw  (q)=-  Ln[det(R„ )]-  tr(R;'R„ ) ,  (2-29) 

then  Equation  (2-28)  becomes 

Lq  (Rn  >  Rs  )=  L^  (q)  -  Ln[det(Rf, )].  (2-30) 

Note  also  that  the  rank  of  the  matrix  is  given  by 

rank(s^)  =  rank(S),withR„  =S„  (2-31) 


Thus,  the  maximization  of  (Rn  ,Rs )  is  the  same  as  the  maximization  of  (q) ,  relative  to 

the  source’s  subspace.  This  question  was  raised  before  and  the  solution  is  known  to  be  given 
by 

max  (q)  =  -  (m  -  q)Ln[ar„  (q)]  +  (m  -  q)Ln[ge„  (q)]  -  Ln[det(R„  )]  C, ,  (2-32) 

where  ar^(q)andge^(q)are  the  arithmetic  and  geometric  means  of  the  (m-q)  lowest 
eigenvalues  of  R^ ,  respectively  and  Ci  is  a  constant  given  by  C,  =  -mLn(7u) .  Then,  using 
Equation  (2-30),  it  can  easily  be  seen  that 

max  Lq  (Rn  ,  Rs  )  =  -  Ln[det(RN )]  -  (m-  q)Ln[ar„  (q)]  +  (m-  q)Ln[ge„  (q)]  -  Ln[det(R„  )]  +  C, 

(2-33) 

However,  recall  from  Equation  (2-27)  that 
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-  Ln[det(R„ )]  =-  Ln[det(R^)]  +  Ln[det(R^)] ,  (2-34) 

which  implies  that  Equation  (2-33)  becomes 

max  Lq  (Rn  ,  Rs )  =  -  (m  -  q)Ln[ar^  (q)]  -F  (m  -  q)Ln[ge„  (q)]  -  Ln[det(R  J]  +  C, .  (2-35) 

Note  in  this  last  equation  that  the  last  two  terms  do  not  depend  on  R^^  or  R^ .  Therefore  in 
the  remainder  of  this  report,  we  will  only  consider  the  following  expression 

Lq  (Rn  )= -  (m  -  q)Ln[ar„  (q)]  -1-  (m  -  q)Ln[ge„  (q)] .  (2-36) 

The  problem  now  is  how  to  maximize  Lq(RN)  (for  a  given  number  of  sources  q), 
relative  to  parameters  defining  the  noise  matrix  R^  .  This  of  course  will  depend  on  the  noise 
parameterization,  which  is  discussed  next. 

2.3  -  Noise  Parameterization 

Consider  first  the  case  of  an  AutoRegressive  (AR)  modeling  of  the  noise.  This  means 
that  noise  received  on  a  sensor  can  be  “predicted”  from  noise  received  on  other  sensors.  This 
assumption  applies  for  the  general  case  that  we  are  considering,  since  usually  we  are  dealing 
with  an  array  of  sensors.  We  will  show  later,  that  actually  all  noise  may  be  described  by 

AR(MA)  modeling.  Assuming  an  AR  type  modeling,  it  can  be  shown  (Gohberg)ll  that  the 
inverse  of  the  covariance  matrix  is  given  by 

r^=:^(a,a;-a,.a;),  (m?) 

where  A,  and  A3  are  two  triangular  (mxm)  Toeplitz  matrices  given  by 


'  *  T.  Kailath  et.  al.,  “Inverse  of  Toeplitz  Operators,  Innovations  and  Orthogonal  Polynomials,”  SIAM  Rev., 
Vol.  20,  pp.  106-119,  Jan.  1978. 
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(2-38) 


■  1 

0 

a,  0 

ap  0 

i 

and  A3  = 

:  Q 

a2 

q  ap  •••  ai  1_ 

a,  3.2  •••  ap 

Note  that  this  formula  is  vahd  for  any  real  stationary  AR  noise  of  order  p  with  coefficients 
a;  ;i  =  l,2,-",p  and  is  the  input  driving  noise  power.  Define  the  matrix  Z’ ,  where 


(j-k)  =  i  ;  l<j<m 
else  ;  1  <  k  <  m  ’ 


(2-39) 


Notice  that  this  matrix  has  the  property  that  Z'"^'  =  Z'Z,  where  Z°  is  the  identity  matrix.  Then 
the  matrices  Aj  and  A3  can  be  written  as 


A,  and  A3  , 


(2-40) 


i=0 


i=0 


where  ag  =  1 ,  by  convention. 

In  the  case  of  an  MA  process,  unfortunately,  there  is  no  exphcit  formulation  of  the 
inverse  of  the  noise  covariance  matrix.  However  the  following  expression  may  be  used 


Rn  = 


ib-X  Xb,y, 

V  1=0  A  >=0  y 


(2-41) 


where  t  denotes  transpose  and  the  bj  ;i  =  l,2,---,p  are  the  MA  coefficients  and  {  }  are 

rectangular  (mx2m)  matrices  defined  by 


Y,a.k)= 


1  if  (k-j)  =  i 
0  else 


(2-42) 


By  combining  Equations  (2-38)  and  (2-41),  we  obtain  an  expression  for  the  parameterization 
of  an  ARMA  process 


(  P  ^ 

(  P  ^ 

II 

z 

Oi 

Ib,Y, 

l^AR 

Sb,Y, 

(i=o  2 

(i=o  y 

(2-43) 
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where  is  the  extended  (2mx2m)  matrix  of  the  AR  process  as  defined  in  Equation  (2- 
37),  but  in  dimension  2m. 

The  noise  model  used  previously  assumes  that  the  noise  is  equal  in  each  sensor  and 
moreover  that  it  has  equal  power.  In  the  case  where  this  assumption  is  not  vahd,  the 
following  noise  model  is  better  suited 

X„+a,X„.,+a,X„.,=|3„w.,  (2-44) 

where  w^  is  a  sequence  of  independent  white  noise  and  is  an  additional  parameter  which 
has  to  be  optimized. 

There  are  several  types  of  parameterizations  for  the  noise  covariance  matrices, 
however,  the  ones  which  we  consider  have  the  advantage  of  defining  a  parameterization  of 
Rf^  (or  R;;,’  )  with  a  small  number  of  parameters. 

In  the  next  section  we  discuss  ways  of  maximizing  the  likelihood  functional  in  the  case 
of  an  AR  model. 

2.4  -  Maximization  of  the  likelihood  functional  in  the  case  of  a  real  AR  model 

The  maximization  of  requires,  in  general,  an  iterative  technique.  A  general  form  of 
the  gradient  algorithm  is  given  by 

Ak+i  “Ak  ~Pk^k ’  (2-45) 

where  Ak=[(Jk  >  »  ■"  »  ap]  is  the  estimated  vector  of  parameters  at  the  k-th 

iteration,  is  the  gradient  vector  at  the  k-th  iteration  and  is  the  step  size  of  the 
algorithm. 

Recall  that  if  Aj'andAj  are  given,  the  inverse  of  the  noise  covariance  matrix  is 
obtained  by 

RN,k=^[A)  .(Af)'  -A3^(A3^)'],whereAf  =Xafz‘  .  (2-46) 

Ct  i=o 
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The  problem  now  is  reduced  to  the  computation  of  the  gradient  vector  G,^  whose  i-th 
component  is  given  by 


G>(i)=— L(q,k), 


(2-47) 


where  L(q,k)  is  the  log-likelihood  functional  conditional  to  q  sources  and  ^ . 

It  can  be  shown  that  computation  of  amounts  to  computation  of  partial  derivatives 
of  the  eigenvalues  of  .  We  will  first  derive  some  properties  in  the  case  of  simple 
eigenvalues  and  compute  the  desired  partial  derivatives. 

•  First  recall  that  R,.,  and  R^  can  be  written  as  R^'^  andR  =  L^R^L”  .  Also  the 

identity  matrix  I  can  be  written  as  .  Therefore 


det(R-;,,R.  -7.,I„)  =  del(L;"L;'L,R.L"  -X„L-"L^)=dM[L;"(R. 
=  det(L;")det(R.,,  -).„I„)det(L;)  =  det(R.,, 


.  (2-48) 


^  1 

This  clearly  shows  that  the  eigenvalues  of  R„  and  Rn  ^R^  are  identical. 


•  Second,  note  that  since  R,^  is  a  positive  matrix,  it  can  be  decomposed  in  triangular 
factors;  i.e.,  R  =  TT”  .  Consider  then  the  same  determinant  as  above 

det(R-;^R,  -k„I„)  =  det(T'"T“R^ JT"  -X„T,-"T")  =  det[T;"(R-„:j 
=  det(T,-")det(R-'^  ->.„I„)de((T^")  =  det(R„\ 

which  means  that  the  eigenvalues  of  R^„  and  T  Rn.rT  are  identical. 

From  the  two  properties  seen  above,  it  can  be  deduced  that  the  computation  of  the 
partial  derivatives  reduces  to  the  computation  of  the  partial  derivatives  (diXj/dajj, 

where  the  Pj's  are  the  eigenvalues  of  the  Hermitian  matrix  T“Rf^\T.  This  result  is  very 

important  in  the  sense  that  the  inverse  of  the  noise  covariance  matrix  is  directly  related  to  the 
AR  coefficients. 


,(2-49) 
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A  classical  result  in  linear  algebra  states  that  given  a  matrix  A  with  its  simple 
eigenvalues  and  its  corresponding  eigenvectors  V^,  the  partial  derivative  (dX^/da-j  is 

given  by 


dX 


da 


J  =:V”  — 

da, 


AV  j ,  where 


da. 


(k,l)=— [A(k,l)]. 


(2-50) 


We  are  now  ready  to  present  all  the  steps  involved  in  the  determination  of  the  maximum 
likelihood  functional. 
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Algorithm 


Computation  of  partial  derivatives  (0  <  i  <  p)  of  N  '  relative  to  the  parameters  {aj } , 

3  1 

A*:  =  —  R-^\  =  —[z\A^y+A\(Z'y-  Z™"  (A3' ) ‘  -  A3'  (Z™-* ) ‘ ] for  1  <  i  <  p .  (2-51) 


‘  da: 


Note  that 


a 


Ai=-2(G,r'R-y,. 


(2-52) 


Computation  of  the  derivatives  of  the  matrices  (A' )  ;  i.e., 


(A^)'  =  —(!'' R'‘,T)  =  T”A^T;0<i<p. 


Computation  of  partial  derivatives  of  the  (simple)  eigenvalues  of  ^ ;  i.e., 


(2-53) 


aa, 


■?^j=(u;)»(A^)'(yf), 


(2-54) 


where  U'  is  the  eigenvector  associated  to  the  eigenvalue  X-'  of  the  matrix  T”R”’^T. 
Computation  of  the  gradient  vector  G,^ ,  defined  by  its  component  G,^  (i) ;  i.e., 


...  ^  m  ^ 

G,(i)=  I  (V,)-^V,-[ar‘(q)]- 


j=q  +  I 


(2-55) 


where 


1  m 

ar'(q)  =  — S^'3. 


(2-56) 


In  the  above  algorithm,  we  have  shown  all  the  steps  necessary  for  the  computation  of 
the  gradient  vector  G^ ,  in  the  case  of  a  real  AR  model.  However,  as  noted  previously,  the 

algorithm  requires  the  knowledge  of  the  eigensystem  of  T^Rj;,’^T.  The  eigenvalues  and 
eigenvectors  of  this  matrix  can  be  computed  exactly  using  standard  algorithms  or  they  may  be 
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estimated  as  discussed  later.  It  is  important  to  note  that  this  algorithm  depends  on  a 
computation  of  a  satisfying  step  size  p  for  the  gradient  algorithm  to  ensure  convergence.  This 
will  be  achieved  by  use  of  a  first  (or  higher)-order  approximation  to  the  change  in  eigenvalues 
of  the  whitened  matrix.  More  precisely,  let  be  the  vector  of  partial  derivatives  of  A,j ;  i.e., 


(A-)'  = 


’da. 


(2-57) 


Then,  to  a  first  order  approximation,  the  j-th  eigenvalue  of  the  whitened  matrix  R„  ^  or 
[RN,k{p)]  'r  with 


RN,k  =:^fAKp)[AKp)r  -A3^(p)[a3^(p)]‘ 


(2-58) 


and 

Ar(p)  =  i[a,-pG.(i)]Z'  (2-59) 

i=0 

is  given  by 

?k';(p)  =  ?t';-pG‘Aj.  (2-60) 

Making  these  substitutions  into  the  likelihood  fiinctional  and  by  means  of  a  uni-dimensional 
method,  an  approximation  for  p  is  determined. 

A  second  order  approximation  may  also  be  used.  In  practice,  however,  a  first  order 
approximation  seems  to  be  sufficient  to  ensure  convergence  of  the  algorithm.  Moreover,  it  is 
less  computationally  extensive. 

2.5  -  Extensions 

In  this  section,  we  extend  the  above  analysis  to  include  different  parameters  and  check 
how  the  technique  performs.  First,  we  extend  the  technique  to  include  what  is  referred  to  as 
Reflection  Coefficients. 
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2.5.1  -  Parameterization  using  reflection  coefficients 


The  parameterization  of  the  matrix  may  also  be  expressed  in  terms  of  reflection 
coefficients.  In  this  way,  the  stability  of  the  AR  model  is  ensured.  Specifically,  the  Levison 
recursion  can  be  written  in  matrix  form  as 


where  A  is  the  [(p+2)x(p+l)]  matrix  defined  as 


A  = 


1  •••  0 

0  •••  1 

0  •••  0 


(2-61) 


(2-62) 


J  is  the  reflection  matrix  ( =1)  and  Ap  is  the  vector  of  AR  coefficients  of  the  noise  model; 
i.e.,  Ap  =[l,a],-",ap].  It  is  then  possible  to  express  Ap^,  in  terms  of  the  reflection 
coefficients  [l,k,,---,kp]  which  leads  to 

Ap«  =(l-k,j)A(l-k,.,j)A.-.(l-k,j)A.l.  (2-63) 

Hence,  the  partial  derivatives  of  the  log-likelihood  functional  can  be  expressed  with  respect  to 
the  reflection  coefficient  as 


(2-64) 


Therefore,  the  partial  derivatives  (3Lq/3kj)  can  easily  be  obtained  since  the  partial 
derivatives  (dL^/da^)  have  been  already  derived  in  the  previous  section  and  (daj/dkj)  can 

be  obtained  from  Equation  (2-63).  The  test  of  stability  for  the  polynomial  A(z)  reduces  to  the 
condition 


2.5.2  -  ARMA  Case 


The  fact  that  an  MA  model  is  included  now  leads  us  to  replace  the  matrix  by  the 
matrix  defined  by 


f  P  ^ 

t 

f  p  ^ 

Sf=-(RN,k)“’- 

Y;  • 

SbfX, 

SbM, 

y; 

U=o  J 

^'=0  J 

(R 


N,k 


(2-66) 


The  ARMA  case  is  obtained  by  derivation  of  Equation  (2-43).  Note  that  the  MA  model  is 
assumed  to  be  minimum  phase. 

2.5.3  -  Multiple  Eigenvalues 


In  the  case  where  some  eigenvalues  have  multiplicity  greater  than  1,  the  formula  giving 
the  eigenvalues’  partial  derivatives  no  longer  holds.  In  this  case,  up  to  the  first  order,  the 
eigenvalues  of  the  perturbed  (R+5R)  are  the  same  as  the  eigenvalues  of  the  matrix 


A=U”  •(R  +  5R)  U, 


(2-67) 


where  U  is  the  matrix  of  an  eigenvector  basis  associated  with  the  multiple  eigenvalue.  This 
leads  us  to  replace  Equation  (2-43)  by  the  computation  of  the  eigenvalues  of  the  matrix 

U“  •  (a*!  )  ■  U .  The  other  steps  of  the  algorithm  remain  the  same. 

2.5.4  -  Complex  Case 

The  fact  that  the  noise  field  is  non  symmetric  with  respect  to  the  array  broadside 
results  in  complex  AR(MA)  coefficients  for  the  noise  model.  However,  nothing  changes  as 
far  as  the  procedure  is  concerned. 

2.6  -  Convergence  Analysis 

In  this  section,  we  study  the  convergence  of  the  gradient  algorithm  defined  in  the 
previous  section.  For  instance,  consider  a  second  order  AR  model.  Its  covariance  matrk  can 
be  computed  using  Equations  (2-37)  and  (2-38).  It  is  then  possible  to  determine  the 
functional  (q)  as  defined  in  Equation  (2-36),  where  ’s  are  the  eigenvalues  of  the 

matrix  Rz!  ^  -R.,  and  ai  and  a2  are  the  AR  parameters.  Note  that  L„  „  (q)  is  the  exact 
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log-likelihood  functional  subject  to  the  condition  on  {apaj.  A  3-dimensional  plot  of  - 
^  (q)  as  a  function  of  the  AR  parameters;  i.e.  aj  anda^,  shows  that  its  maximum  is 

attained  at  the  exact  values  of  the  AR  parameters  (see  Figures  1-1  and  1-2).  However,  notice 
that  the  function  is  not  concave,  however  the  following  properties  do  apply  to  the  function 

La,.a,(q): 

1.  The  log-likelihood  function  is  negative.  It  equals  0  if  and  only  if  the  noise  model  is 
estimated  perfectly. 

2.  The  gradient  vector  of  the  log-likelihood  function  is  null  if  and  only  if  A  is  equal  to 
Aq  ,  which  represents  the  exact  values  of  the  coefficients. 

In  previous  sections,  we  assumed  that  the  number  of  sources  (q)  and  the  AR  order  (p) 
are  fixed.  Next,  we  will  look  at  the  consequences  of  misadjustment  of  these  parameters. 

2.6.1  -  Number  of  sources  (q) 

Using  property  1  above,  it  is  possible  to  overestimate  the  number  of  sources  q,  as  long 
as  the  order  of  the  AR  model  p  in  less  than  (m-q)  without  degradation  in  the  asymptotic 
case.  In  practice,  the  covariance  matrix  is  estimated,  therefore  the  overestimation  of  q  will 
lead  to  shghtly  inferior  performances  of  the  method.  However,  the  degradation  is  quite 
acceptable.  The  question,  however,  is:  what  would  be  a  good  strategy  for  choosing  the  value 
of  q.  Recall  that  the  classical  information  criteria  such  as  those  developed  by  Akaike  and 
Rissanen  do  not  provide  a  satisfying  estimation  of  the  number  of  sources.  This  is  due  to  the 
fact  that  they  make  use  of  the  eigenvalues  of  the  estimated  covarianee  matrix  R,^  and  they  do 
not  separate  between  the  uncorrelated  sources  and  (highly)  correlated  sourees.  To  solve  these 
problems,  we  present  the  following  solutions 

a)  Source  over-determination: 

Trying  conseeutive  values  of  q  and  then  choose  the  value  which  maximizes  the  log- 
likelihood  function  will  lead  us  to  the  right  answer. 
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b)  Modifications  of  information  criteria: 

The  information  theoretic  criteria  are  based  upon  log-likelihood  determination  and  its 
statistical  behavior.  These  statistics  are  characterized  by  a  number  of  “free”  parameters. 
The  fact  that  we  have  an  AR  noise  model  enables  us  to  add  some  “free”  parameters. 
More  specifically,  assume  that  Rissanen’s  criterion  is  used,  in  which  case,  the  following 
function  is  considered 

MDL(p,  q)  =  -L(p,q)  •  [p{2m-  p)  +  q  +  2] ,  (2-68) 

where  N  is  the  number  of  independent  samples,  p  is  the  number  of  sources  and  q  is  the 
noise  model  order.  L(p,q)  is  the  likelihood  function  conditional  on  the  following' 
hypotheses:  p  sources  and  noise  model  order  q. 

c)  Statistical  properties  of  the  eigenvalues  of  a  “corner”  of  matrix  R  (MA  case): 

This  original  approach  has  been  presented  by  Fuchs  ^  2  and  gives  an  accurate  estimate  of 
the  number  of  sources  in  the  presence  of  correlated  noise. 

d)  Use  of  a  State-Space  Approach  for  Sensor  Outputs  Modeling  (ARMA  Case): 

Using  an  information  criterion  such  as  the  predictive  efficiency  criterionl^,  it  is  possible  to 
obtain  satisfactory  estimates  of  the  number  of  sources  without  a  priori  knowledge  about 
the  noise  model. 

2.6.2  -  Noise  Order  Model 

In  the  asymptotic  case,  noise  model  over-determination  leads  to  very  slight 
degradation  of  the  results  in  terms  of  noise  spatial  density.  We  will  show  later,  through 
computer  simulations  that  it  is  possible  to  overestimate  the  order  of  the  noise  model,  without 


12  J.  J.  Fuchs,  “Estimation  du  Nombre  de  Sinusoids  Dans  du  Bruit  Colore,”  Actes  du  Onzieme  Collogue 
GRETSI.  Vol.  1,  pp.  197-200,  June  5,  1987. 

12  K.  S.  Arun  and  S.  Y.  Kung,  “Generalized  Principal  Components  analysis  and  its  applications  in 
Approximate  Stochastic  Realizations,”  in  Modeling  and  Applications  of  Stochastic  Processes,  U.  B.  Desai,  Ed. 
Boston:  Kluwer  Academic,  1986,  pp.  75-105. 
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dramatic  effects.  As  shown  previously,  it  is  obvious  to  note  that  the  parameters  p  and  q  may 
be  estimated  simultaneously.  However,  a  separate  estimation  of  p  seems  rather  difficult. 

Another  point  to  consider  is  the  stability  of  the  AR  model  obtained  by  the  maximizing 
of  the  log-likelihood  function.  It  is  possible  to  compute  the  roots  of  A^(z) ,  where  A|^(z)  is 
the  polynomial  associated  with  A^ ,  and  to  adjust  the  step  size  of  the  gradient  method. 
However,  it  can  be  shown  that  the  parameterization  by  reflection  coefficients  is  best  suited  in 
this  case.  It  is  interesting  to  note  that  in  practice,  if  the  noise  poles  approach  the  unit  circle, 
then  the  noise  covariance  matrix  tends  towards  singularity.  This  leads  to  the  log  likelihood 
functional  approaching  infinity.  However,  if  the  step  is  not  very  important,  this  avoids  any 
stability  problems. 

2.6.3  -  Estimation  of  Source  Parameters 

Assume  that  an  estimation  of  the  AR  parameters  has  been  obtained  after  several  runs 
of  the  algorithm  described  earlier.  Let  A‘ =|l  ,  a,  ,  •••  ,  Sp}  be  this  estimate. 

An  apphcation  where  the  above  problem  arises  would  be  the  estimation  of  Angles  Of 
Arrival  (AOA)  of  some  sources  as  well  as  their  powers  at  a  given  frequency.  The  derivation 
of  source  powers  estimation  is  a  fundamental  tool  for  the  minimization  of  the  whiteness 

functional.  The  whitened  matrix  R„  given  by  R„=L-R-L”,  with  R^  =L-L”,  is 
considered  for  source  bearing  estimation.  Note  that  R^,  is  not  generally  a  Toephtz  matrix 
(even  if  R,^  is  Toeplitz),  however,  we  can  remedy  this  problem  through  a  filtering  process. 
The  steering  vectors  Dg ,  corresponding  to  bearing  0,  becomes  whitened;  i.e.,  (L.  Dg). 

The  MUSIC  method  for  AOA  estimations  consists  of  computing  the  sine  of  the 
steering  vector  (L.  Dg)  on  the  source  subspace  (space  spanned  by  the  eigenvectors 

corresponding  to  the  largest  eigenvalues  of  R^ , 


1 

sin'(e) 


||L-Deir 

|(i-n)LD, 


|L-De 

(nq.LD, 


(2-69) 
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where  ||  ||  is  the  EucUdean  norm,  11  is  the  projection  matrix,  11  =  V  •  V“  and  v(mxq)  is  the 
matrix  constituted  by  the  greatest  eigenvectors  of  R„ . 

Pisarenko’s  method  has  also  been  considered.  Let  jPj  ,  be  the 

estimated  angles  of  the  sources.  Consider  the  (mxm)  associated  matrices;  i.e.. 


where  a j  ^  is  the  1-th  row,  k-th  column  element  of  the  matrix  |d§  •  D?  |  for  i  =  1, 2,  •  •  • ,  s  and 


1  =  l,2,-",(m-p) .  All  the  sub-matrices  are  Toephtz,  but  non-Hermitian.  Then  the 
transformed  covariances  of  sources  are  given  by 

s.(e,,,)=A".s(e„)-A-s4e„i).  (2-71) 

The  covariance  matrices  of  the  sources  s„(ej,lj  correspond  to  the  covariances  of  the 

sources  after  whitening  by  the  inverse  filter  A  as  will  be  seen  later.  Covariances  of 
transformed  array  outputs  (1)  are  defined  by  the  same  method  as  previously  shown,  using  a 
Toephtz  estimate  of  the  covariance  matrix  of  the  outputs.  This  estimate  is  obtained  by 
averaging  along  the  diagonal  and  is  the  same  as  an  orthogonal  projection  on  the  Toephtz 
subspace.  Pisarenko’s  method  can  then  be  used  for  estimating  the  power  of  the  sources. 
Assuming  spatial  whiteness  of  the  additive  noise,  the  estimates  of  the  source  powers 

|g,  •  sgq  }  are  the  solutions  of  the  following  over-determined  linear  system 


(2-72) 
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It  is  important  to  observe  that  the  procedure  described  above  is  similar  to  the  case  of 
uncorrelated  sources  which  means  that  the  source  matrix  is  given  by 


S  =  ig,D.^De%  (2-73) 


and  theoretically  is  independent  of  the  noise  level  since  it  uses  covariances  (1)  with  1  >  1 . 

2.7  -  Plane  Wave  Hypothesis,  Whiteness  Functional 

In  this  section,  we  will  focus  on  the  case  of  a  linear  equispaced  array  of  sensors.  We 
will  derive  the  expressions  for  the  likelihood  functional  in  the  case  of  an  AR  system,  using 
Pisarenko’s  technique  for  harmonic  retrieval. 

We  assume  that  the  covariance  matrix  R  is  Toephtz  and  that  the  source  parameters  are 
perfectly  known  and  define  the  function 

J  =  X|r(l)-r,(l)f  ;(L>1)  .  (2-74) 

1=1 


Note  that  J  is  null  in  the  case  of  white  noise.  In  the  case  where  the  noise  is  correlated, 
it  is  possible  to  use  whitening  to  reach  a  similar  conclusion. 

We  will  only  consider  the  case  of  an  AR  model.  Let  A  be  the  vector  of  coefficients  of 
this  model  and  let  [^A(z“')j  be  the  associated  filter.  The  whitened  filter  is  A(z“’)  and  the 
covariances  of  whitened  data  can  be  written  as 


V.i=0 


Vi=0 


k  where  a„  =  1. 


(2-75) 


Therefore, 

rw(l)  =  A”-R,-A,  (2-76) 


where  R,  is  the  (p-i-l)x(p+l)  matrix  given  by 
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r(l-l) 

r(l) 


(2-77) 


r(l) 

rd  +  l) 


-  r(l-p) 
r(l-p-l-l) 


[r(l-l-p)  r(l-l-p-l)  •••  r(l)  J 
Note  that  the  above  matrix  is  Toeplitz  but  non-Hermitian.  Define  the  function 

J.(A)  =  t|r.(l)-r,(l)f  ,  (2-78) 

1=1 


where  r^(l)  was  defined  earlier  and  rj(l)  are  the  exact  covariances  of  sources  after  the 
mapping  A.  Then  when  A  is  equal  to  A^ ,  where  A^  is  the  exact  model,  we  can  show  that 


r.(l)  =  r.(l)  +  E 


Anfz  ')  Ao(z  ’) 

Ao(z  )  Ao(z  ) 


=  r,(l), 


where  e(t)  is  the  input  white  noise  sequence  for  the  AR  model.  Therefore 

I.(Ao)  =  0- 

For  practical  applications,  one  uses  the  whiteness  function  defined  by 

Jw(A)  =  IK(l)-f,(l)|\ 


(2-79) 


(2-80) 


(2-81) 


1=1 


where  f^(l)  was  defined  earher,  but  using  the  available  covariance  matrix  R,  instead  of  the 
exact  matrix  R, ,  and  f^(l)  is  estimated  using  Pisarenko’s  method  or  any  other  high  resolution 
technique  applied  to  the  whitened  data  in  the  following  manner 


m 

fs(l)  =  exp{-  }’ 

i=l 

where  m  is  the  number  of  sensors,  gj  and  f^  are  the  power  and  spatial  frequency  of  the  i-th 
source,  all  of  which  were  estimated  using  a  high  resolution  technique.  Note  that  L  is  usually 
chosen  to  be  L  =  (m  -  p  - 1) . 


Ill 


We  will  now  deal  with  the  minimization  of  the  fimction  defined  by  Equations  (2-81) 
and  (2-82).  This  minimization  is  done  in  terms  of  AR  coefficients  vector  A  only,  but  is 
relative  to  the  terms  (1)  and  (1) .  The  general  procedure  is  summarized  below; 

•  Iteration  0,  (Starting  Parameters): 

A‘  =[l,0,-,0] 

g°  and  f  °  obtained  by  spatial  analysis  of  R 

rw(l)  =  f(l) 

•  Iteration  k: 

— k  “ [^’^1 

gf  and  fj'‘  obtained  by  spatial  analysis  of  R^ 
f,^(l)  =  A«-R,-A. 

Note  that  the  number  of  sources  may  be  corrected  at  each  iteration. 

•  Iteration  (k+l):  Computation  of  the  gradient  vector  Gj  ) 

=  A^-P^Gj(a^),  where  is  the  gradient  step  size . 

2.7.1  -  Whitenin2  Invariance  Properties 

Recall  that  r^  (1)  can  be  expressed  in  matrix  form  as  R^  =  A  •  R^  •  A “  ,  where 

1  aj  ap  0  . 

0  0 
i  0 

1  0  1  a,  •••  ap  0  ■■■ 

Here  A  is  the  directional  matrix  and  therefore  should  not  be  confused  with  the  gradient  vector 
discussed  in  section  2.7.  Note  also  that  the  matrix  R^  defined  earlier  is  a  Toeplitz  positive 


0 


0 


(2-83) 
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matrix.  Moreover,  we  have  R,,(l,l)  =  This  matrix  transformation  aUows  us  to  study 

the  effect  of  whitening  on  the  covariance  matrix. 

Next,  we  will  study  the  effect  of  whitening  on  a  plane  wave  whose  associated  steering 


Dl  =  1  ,  ej“  ,  •••  ,  e 


vector  is  given  by 


It  can  then  easily  be  seen  that  the  covariance  matrix  is  given  by 

R.(0)  =  ADeD»A«. 


(2-84) 


(2-85) 


Note  that  the  product  ADg  is  given  by 


ADe  = 


l  +  a,eJ“-h...-t-apeJ''’-''“ 
e^“ -Ha,e'j“-l-...-f-apej'’“ 


which  can  also  be  written  as 


ADe  =  (l  +  a,e -f-. .  .-t-ape^^"''’"  )£e”  * 


(2-86) 


(2-87) 


where  Dg  is  and  (Lxl)  vector  obtained  by  selecting  the  first  L  components  of  the  vector  Dg . 
Hence,  it  can  easily  be  shown  that  the  covariance  matrix  can  be  expressed  as 

R.(e)  =  q(e)D,.D;,  (2-88) 


where 


q(0)=  l-t-a,ej“4-...-i-ape 


j(p-l)a| 


(2-89) 


The  above  properties  lead  to  the  following; 

•  Let  R^be  a  Toephtz  covariance  matrix,  then  R^=A  R^  A”  is  a  Toephtz 
matrix. 
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•  The  procedure  transforms  a  steering  vector  into  another  steering  vector 
corresponding  to  the  same  angle  of  arrival. 

In  practice,  the  exact  angle  of  arrival  of  the  sources  are  unknown.  The  above 
properties  state  that  only  these  angles  are  invariant  under  the  whitening  procedure.  These 
angles  may  be  estimated  using  any  high  resolution  technique.  It  is  weU  known,  however,  that 
these  estimates  will  never  be  exact.  Using  the  whiteness  fimction,  the  intent  is  to  define  an 
iterative  method  which  will  minimize  the  differences  between  the  off-diagonal  terms  of  the 
whitened  matrix  and  the  corresponding  terms  of  the  source  matrix. 

2.7.2  -  Calculation  of  the  Gradient  Vector  of  the  Whiteness  function 

Using  Equations  (2-81)  and  (2-83),  we  can  show  that  the  gradient  vector  is  given  by 

G,(A)  =  2Re|x[R,A-ar.(l)][A"  R,"  ^-^(1)]!,  (2-90) 

where  ^^(l)  is  the  gradient  vector  of  r3(l,A)  relative  to  A.  The  hardest  part  of  the  problem 
is  to  determine  this  last  quantity.  Note  that  Gr^.  (1)  is  defined  by 

Qrs(l)  =  Za(A)'e-«  -jJtl  XgiS-W  ■Gf,(A).  (2-91) 

i=l  i=l 

Note  that  the  computation  of  Grj(l)  amounts  to  the  computation  of  Gfi(A)  and 
Ggj  (a)  .  This  will  be  achieved  in  two  steps.  First,  we  compute  the  partial  derivatives  of  the 
eigenvectors  of  R^.  Then,  derivatives  of  the  source  parameters  are  computed  using  a 
perturbation  analysis  of  a  high  resolution  technique. 

1.  Partial  Derivatives  of  Eigenvectors  of  R„ 

Denote  by  R^  the  (L+l)x(L-i-l)  Toeplitz  matrix  defined  by 

RA(i-j)  =  fA(i-j)>  (2-92) 

with  fA(l)  =  A” -R,  •  A  and  {yj,-",UL+,;A,,  is  the  eigensystem  of  R^  .  In  the 

case  of  a  simple  eigenvalue,  it  can  be  shown  that 
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U-k  -  2-i 


aa, 


j=l 

j#k 


A-k  Aj 


Uj;  (n  +  l<k<L  +  l.  (2-93) 


The  matrix  (aR  A /aaj )  is  a  Toeplitz  matrix  which  can  be  defined  by  its  first  row  as 


^R*(m,l)  =  E^R„.,  -  A  +  A"  .R„_,  E, , 

aa, 


(2-94) 


with 


iKi) 


ri;l  =  i 

[0;else 


(2-95) 


2.  Derivatives  of  Source  Bearings 

Recall  that  the  i-th  source  is  characterized  by  its  power  and  spatial  frequency  gj  and  f j , 
respectively.  The  problem  now  consists  on  computing  their  partial  derivatives  with  respect  to 
the  parameters  {aj}.  In  the  case  where  the  MUSIC  algorithm  is  used,  recall  that  the 
estimates  of  the  angles  of  arrival  are  obtained  by  minimizing  the  projection  of  vector  Dg  on 
the  noise  subspace:  n  =  U-U”,  where  U  =  [U„^,  ,  •••  ,  U^]  (“lowest”  eigenvectors  of 

Rw). 

The  projection  of  Dg  on  the  noise  sub-space  is  given  by: 

n(e)  =  Dg»nDg.  (2-96) 


Letz‘=[l  ,  Zj  ,  ,  z^^jwithZj  =  exp|-j7tfi  j.  Then  the  derivative  of  this  vector 

with  respect  to  the  spatial  frequencies  is  given  by 

-iz;=-j7i[0  ,  z,  ,  -  ,  LzH.  (2-97) 

Then  an  expression  of  the  partial  derivatives  of  the  spatial  frequencies  with  respect  to  the  AR 
coefficients  is  expressed  as 
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d  ah 


^  H 


H  3 


3a, 


z  n-Zj+z  -n— z.+z  •— z.  =o 


3a, 


So  far  we  have  shown  how  to  compute  the  partial  derivatives  of  the  eigenvectors  using 
equation  (2-93),  the  spatial  frequencies  fj  using  a  high  resolution  technique  and  the  partial 

derivatives  ^3fi/3aij.  If  the  projection  of  the  whole  noise  subspace  is  considered,  the 
derivative  of  the  projector  n  is  given  by 


Re 


3  -  _  1 


H 

T"  * 

■  z, 

•  z,  -u 

[3a,  J 

-J 

3a, 


K 


Im[(u"Z,)(z'u)} 


(2-98) 


and 


3a, 


■n(9)= 


v3a, 


J 


n-D'e+D” 


3a| 


n'e+D'e-n 


J 


v3a, 


(2-99) 


The  next  step  will  involve  computing  the  partial  derivatives  of  the  projector  11.  To  do  this, 
we  take  a  second  order  of  IT  as  given  by 


Hg  =(Uo  +6U  +  5^U)-(l-||6u|f)-(u”  -hbU”  +5^U”), 


(2-100) 


where 


u,s(u,),+(8y,)+(yu,)  (2-101) 

and 

(Uo-f5U  +  5'u)”  •(Uo+5U-i-6'u)  =  I  +  ||6uf .  (2-102) 

From  Equations  (2-100),  (2-101)  and  (2-102),  we  can  get  a  first  order  expansion  of  Ilg  as 

Hg  =n  +  U-5U”-i-5UU“.  (2-103) 

Therefore 
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— n  =  u-  -U",  (2-104) 

aa,  l^da,  j  {^a^ 

from  which  we  obtain 


It  was  possible  to  have  used  a  second  order  expansion.  However,  in  most  practical  situations, 
a  first  order  expansion  yields  good  results. 

3.  Derivatives  of  Source  Powers 

Let  F  be  an  (Lxq)  matrix,  obtained  from  elementary  theoretical  covariances  of 
estimated  source;  i.e., 

■s(e„i)  -  s(0,,i)l 

F=  i  •••  :  with  s^0;^,lj  = -sin ^Tclfijj.  (2-106) 

s(0„l)—  s(0,,l)J 

Moreover,  let  J  be  the  (Lxq)  vector  obtained  from  the  imaginary  parts  of  whitened  output 
covariances;  i.e., 

r  =Im[fA(l),-",fA(L)],  (2-107) 

where  fA(l)  and  A  were  defined  previously.  Also,  let  F  be  the  vector  of  source  powers 
defined  as 

r'=[g,  ,  -  .  g,].  (2-108) 

In  the  white  noise  ease,  the  vector  F  is  the  solution  to  the  following  linear  system 

FF  =  1-  (2-109) 

The  gradient  vector  G,  where 
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is  obtained  by  differentiating  Equation  (2-109)  and  is  given  by 


fail 

•J-hF~‘  • 

y 

(2-111) 


where 


- F  is  obtained  from  Equations  (2-93)  and  (2-105),  - — J  is  given  by  Equations  (2- 

da,  oa, 

94)  and  (2-95)  and  F'’  is  the  pseudo-inverse  matrbc  of  F.  The  vectors  G-  (A)  are  obtained 

from  the  vectors  G, . 
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Computer  Simulations  and  Results 

In  this  section,  we  present  the  results  of  the  computer  simulations  which  we  carried 
out.  The  goal  of  the  simulations  is  to  demonstrate  the  effectiveness  of  the  proposed 
techniques.  Note  that  at  this  stage  only  the  Parametric  beamforming  algorithm  based  on 
LeCadre's  technique  was  evaluated.  Note  also  that  only  spatial  processing  was  considered. 

Recall  that  the  objective  of  the  proposed  method  is  to  determine  the  number  of  targets 
present  in  the  environment.  Therefore  the  detection  problem  which  we  undertook  in  this 
effort  is  different  from  the  classical  detection  problem  which  we  usually  refer  to  in  the 
hterature  as  a  decision  making  about  the  presence  or  absence  of  a  target.  LeCadre's  technique 
is  beheved  to  improve  the  detection  of  targets,  especially  for  cases  of  low  signal-to-noise 
ratios.  This  is  of  great  interest  to  us  since  in  most  apphcations  the  targets  are  buried  in  the 
environment  consisting  of  clutter,  jammers  and  thermal  noise.  Also,  by  being  parametric,  the 
technique  shows  great  promise  in  the  overall  detection  scheme  such  as  computational 
efficiency. 

3.1  -  Likelihood  Functional 

The  first  goal  of  the  computer  simulations  was  to  duphcate  the  results  presented  in  [1], 
and  then  extend  the  technique  to  more  complicated  scenarios.  Preliminary  computer 
simulated  results  show  that  the  technique  performs  as  expected,  since  we  were  able  to 
duphcate  most  of  the  results  presented  by  J.  P.  LeCadre.  For  instance,  we  showed  that  the 
likelihood  functional  defined  in  the  paper  is  minimized  (or  maximized,  depending  on  which 
quantity  we  are  using)  with  respect  to  the  Auto  Regressive  (AR)  coefficients  of  the  noise 
model  and  this  was  done  for  a  fixed  and  known  number  of  sources.  The  scenario  used  in  this 
case  was  as  follows.  It  consisted  of  three  (q  =  3)  targets  with  parameters  as  shown  in  the 
table  below: 
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Target 

0 

SNR  (dB) 

1 

40 

-5 

2 

50 

-10 

3 

70 

-15 

Table  3-1.  Target  Parameters 


The  AR  model  was  assumed  to  be  of  order  two  (p  =  2)  with  the  following  AR  parameters: 

ao  =  l,  a,  =-0.9,  a2=0.2. 

We  considered  a  linear  uniformly  spaced  array  consisting  of  14  sensors  (  m  =  14),  spaced  at 
half  wavelength  d  =  ?i/2,  such  that  a)d/c=27t(d/>.)=7C.  500  samples  (snapshots)  were 
used  in  the  computation  of  the  sample  covariance  matrix.  The  value  of  the  likelihood 
functional  was  found  to  be  L  =  0.1582.  We  varied  the  values  of  aj  from  -2  to  +2  and  the 
value  of  a2  from  -1  to  +1.  Note  that  throughout  the  analysis  as  well  as  the  computer 
simulations,  the  value  of  ao  is  assumed  to  be  equal  to  unity  (ao  =  1).  Figures  3-1  and  3-2  show 
the  variations  of  the  likelihood  functional  as  a  function  of  ai  and  a2. 
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Likelihood  Functional 


Figure  3-1.  Likelihood  Functional  as  a  Function  of  ai  and  a 


Figure  3-2.  Likelihood  Functional  (in  dB)  as  a  Function  of  ai  and  a? 
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From  the  above  figures,  it  can  be  seen  that  the  minimum  of  the  likelihood  functional  is 
attained  at  the  exact  values  of  the  AR  parameters,  namely  ai  and  ai.  This  shows  that  the 
algorithm  performs  exactly  as  expected. 

3.2  -  AR  Parameters  Estimation 

Another  aspect  of  the  technique  is  to  accurately  estimate  the  order  as  well  as  the 
coefficients  of  the  AR  model  used  to  characterize  the  noise.  In  the  next  simulation,  the  same 
scenario  as  previously  described  is  used.  We  also  assume  that  the  order  of  the  AR  model  and 
the  number  of  sources  have  accurately  been  estimated.  The  following  figures  show  the  AR 
coefficients  estimates  as  a  function  of  the  step  size  (p)  used  and  the  number  of  iterations 
needed  for  convergence.  Note  that  in  all  cases,  we  have  assumed  a  total  number  of  50 
iterations.  However,  as  can  be  seen  from  these  figures,  convergence  was  attained  much  faster 
than  50  iterations  depending  on  the  value  of  the  step  size.  In  the  figures  below,  the  value  of  p 
was  varied  from  0.4  to  0.7  with  increments  of  0.05.  As  can  be  seen  from  these  figures,  the 
bias  errors  seem  to  have  a  sinusoidal  relation  to  the  step  size. 
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Figure  3-3.  AR  Parameter  Estimates  (step  size  =  0.4) 
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From  Figures  (3-3)  to  (3-9),  the  following  observations  were  made.  With  a  step  size 
of  0.4,  convergenee  was  attained  after  5  iterations  (very  fast).  However,  there  is  a  bias  of 
about  0.02  for  both  ai  and  a2.  This  is  not  bad,  depending  on  whether  the  technique  is  robust  to 
errors  in  the  AR  parameters  or  not.  This  will  discussed  later.  Note  that  with  a  step  size  of 
0.45,  convergence  was  attained  after  15  iterations,  however  with  stUl  a  bias  of  about  0.18  for 
ai  and  0.05  for  a2.  With  a  step  size  of  0.5,  convergence  was  attained  after  20  iterations, 
however  again  with  a  bias  of  approximately  0.05  for  both  ai  and  a2.  With  a  step  size  of  0.55, 
convergence  was  attained  after  15  iterations  for  both  ai  and  a2,  with  a  very  large  bias  of  about 
0.25  for  ai  and  a  small  bias  of  about  0.02  for  a2.  With  a  step  size  of  0.6,  convergence  was 
attained  after  15  iterations,  however  with  a  very  small  bias  of  approximately  0.0005  for  aj  and 
a  bias  of  about  0.02  for  a2.  With  as  step  size  of  0.65,  convergence  was  attained  after  5 
iterations  for  a2  with  a  bias  of  0.002  and  after  10  iterations  for  a\  with  a  bias  of  approximately 
0.003.  When  the  value  of  the  step  size  was  0.7,  convergence  was  attained  after  12  iterations 
for  both  ai  and  a2,  however  with  bias  of  about  0.0015  for  ai  and  about  0.25  for  a2.  From  the 
above  discussion,  it  is  clear  that  the  best  results  were  obtained  when  p  was  equal  to  0.6.  This 
value  should  be  the  one  to  be  used  in  further  processing  of  the  received  signals. 

3.3  -  Number  of  Signal  Estimation 

Next,  we  show  that  an  extension  of  the  MDL  algorithm  can  also  be  used  to  determine 
the  number  of  sources  and  the  order  of  the  AR  model.  Computer  simulations  have  shown  that 
the  minimum  of  the  functional  was  attained  at  the  exact  values  of  these  parameters.  Again,  as 
in  the  previous  case,  we  considered  the  same  scenario  as  before.  We  applied  the  modified 
MDL  criterion  as  defined  earlier.  Note  that  these  results  were  achieved  through  whitening  of 
the  covariance  matrix  of  the  received  vector.  Table  2  shows  the  results  of  the  computer 
simulation. 
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Process  Order 

0 

1 

2 

3 

4 

5 

Number  of  Signals 

1 

45.2716 

29.4954 

3.0568 

3.4702 

3.8657 

4.2586 

2 

45.2431 

29.4668 

3.4417 

3.8372 

4.2301 

3 

45.2336 

29.4573 

^^9 

3.4322 

3.8277 

4.2206 

4 

45.3001 

29.5239 

3.0853 

3.4987 

3.8942 

4.2872 

5 

45.4997 

29.7235 

3.2849 

3.6983 

4.0939 

4.4868 

6 

46.5264 

30.0113 

3.6747 

4.0881 

4.4836 

4.8765 

Table  3-1.  Values  of  New  MDL  Criterion 


As  expected,  the  minimum  value  of  the  new  MDL  criterion  was  obtained  when  we 
reached  the  exact  value  of  the  AR  order  model  and  the  exact  number  of  sources.  From  this 
table  we  deduce  both  the  order  of  the  AR  model  (p  =  2)  and  the  number  of  sources  (q  =  3)  for 
subsequent  analysis.  It  is  important  to  note  here  the  efficacy  of  the  algorithm  in  estimating 
these  two  important  parameters.  Without  the  knowledge  of  these  two  parameters,  it  would  be 
extremely  hard  to  process  the  data  any  further,  especially  if  we  are  interested  in  estimating  the 
locations  of  the  sources. 

3.4  -  Angle  of  Arrival  Estimation 

With  the  above  knowledge  in  hand,  we  applied  two  different  algorithms  to  locate  the 
targets  involved  in  the  scenario,  namely  MUSIC^^  and  ESPRITI^  .  It  is  important  to  note 
here  that  two  different  simulations  were  carried  out. 

3.4.1  -  Signal  Plus  Noise  Case  With  Known  Whitening  Functional 

First,  we  assumed  that  the  received  signal  was  of  the  following  form: 

Y=AS+CN, 


'4  R.  O.  Schmidt,  “Multiple  Emitter  Location  and  Signal  Parameter  Estimation,”  IEEE  Trans.  Antennas  and 
Propagation.  Vol.  AP-34,  pp.  276-280,  March  1986. 

*5  A.  Paulraj,  R.  Roy  and  T.  Kailath,  "Estimation  of  signal  parameters  via  rotation  invariance  technique- 
ESPRIT,"  Proc.  19-th  Asilomar  Conf.  on  Circuits.  Systems  and  Computers,  pp.  83-89,  Pacific  Grove,  CA, 
Nov.  1985. 
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where  A  is  the  directional  matrix  (representing  the  array  manifold),  S  is  the  impinging  signal 
vector  and  C  is  the  noise  driving  term  and  N  is  the  additive  white  noise  vector  assumed  to  be 
independent  from  the  signals.  As  described  earlier,  C  was  obtained  from  the  Cholesky 
decomposition  of  the  AR  covariance  matrbc  and  thus  provides  the  specified  cross-channel 
correlation.  We  refer  to  this  matrix  as  the  whitening  functional.  At  first,  we  assumed  this 
matrix  to  be  known  a  priori,  which  means  that  the  AR  covariance  matrix  was  computed 
exactly  knowing  the  exact  values  of  the  AR  parameters. 

3.4.1.1  -  MUSIC  Operator 

Figure  3.10  shows  the  result  of  the  MUSIC  algorithm.  The  solid  line  shows  the 
transformed  MUSIC  algorithm,  while  the  dotted  line  shows  the  un-whitened  MUSIC 
algorithm.  It  is  very  clear  from  the  above  figure  that  whitening  does  improve  the  location  of 
the  signals.  Without  any  whitening,  one  sees  that  MUSIC  was  unable  to  resolve  all  targets. 
Recall  that  the  signals  are  completely  buried  in  the  clutter  (SNR  of -5,  -10  and  -15  dB).  These 
results  clearly  show  the  need  for  the  whitening  functional. 
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3.4.1.2  -  ESPRIT  Operator 

In  the  previous  section,  we  concentrated  on  using  LeCadre’s  technique  in 
conjunction  with  the  MUSIC  algorithm.  It  is  well  known  that  one  of  the  disadvantages  of 
MUSIC  is  the  computational  load  required  to  execute  the  algorithm  as  well  as  the  amount  of 
memory  needed  to  store  the  array  manifold.  For  these  reasons,  we  shifted  our  attention  to 
another  technique  which  is  computationally  less  intensive  and  does  not  require  a  lot  of 
memory  space.  This  technique  is  referred  to  as  ESPRIT  (Estimation  of  Signal  Parameters  via 
a  Rotational  Invariance  Transformation).  In  this  technique,  two  sub-arrays  are  formed.  Then, 
we  form  the  covariance  matrix  of  the  first  sub-array,  call  it  Mi,  and  the  cross-covariance 
matrix  of  sub-array  1  and  sub-array  2,  call  it  M2.  It  can  then  be  easily  shown  that  the 
information  of  the  locations  of  the  targets  is  contained  in  the  rank  reducing  values  of  the 
matrix  pencil  (Mi-^iMa).  These  rank  reducing  values  are  shown  to  be  the  generalized 
eigenvalues  of  the  same  matrix  pencil.  In  the  initial  phase,  we  assumed  that  the  whitening 
transformation  is  known  a  priori.  We  have  selected  the  same  scenario  as  described  previously 
in  which  we  considered  an  array  consisting  of  14  sensors  spaced  at  half  wavelength  and  3 
targets  located  at  40°,  50°  and  70°,  with  an  SNR  of  0  dB,  each.  The  AR  process  was  stiU 
assumed  to  be  of  order  2  with  parameters  1,  -0.9  and  0.2.  The  following  show  the  mean,  the 
variance  and  the  mean-squared  error  of  the  estimates  obtained  after  50  trials. 


0 

0 

II 

cd“ 

02=50” 

II 

0 

Mean 

41.0609 

51.3136 

70.6455 

Variance 

0.0036 

0.0030 

0.0057 

MSE 

56.2737 

86.2731 

20.8348 

Table  3-3.  ESPRIT  Estimates 


We  can  also  visualize  these  results  by  looking  at  the  angle  estimates  as  a  function  of  the 
number  of  iterations  and  the  histograms  of  the  estimates. 
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Angle  =  50  Degrees 


130 


Number  of  Occurences  Angle  -  70  Degrees 
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Figure  3-13.  Angle  =  70°  as  a  function  of  iteration  number 


Angle  =  40  Degrees 


Figure  3-14.  Angle  =  40°  (Histogram) 
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Number  of  Occurences  Number  of  Occurences 
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Even- though  the  angle  estimates  have  a  small  bias,  the  varianees  of  these  estimates  are 
very  small.  This  clearly  shows  the  effectiveness  of  the  proposed  technique. 

3.4.2  -  Signal  Plus  Noise  Case  With  Estimated  Whitening  Functional 

In  this  case,  we  estimate  the  whitening  likelihood  functional  using  the  method  derived 
in  section  2.  Recall  that  the  first  step  in  this  method  is  to  estimate  the  order  of  the  AR  model 
used  to  describe  the  noise  and  the  AR  parameters,  along  with  the  number  of  signals.  Once 
these  parameters  have  been  estimated,  we  generate  the  AR  covariance  matrix  from  which  the 
whitening  functional  is  derived.  This  functional  is  then  applied  to  the  covariance  matrix  of  the 
received  data. 

3.4.2.1  -  MUSIC  Operator 

In  this  section,  we  ran  the  same  scenario  as  described  previously.  However,  in  this 
case  the  whitening  functional  was  estimated  whereas  in  the  previous  case  the  true  whitening 
functional  was  used.  The  following  figure  shows  the  result  of  the  MUSIC  operator.  Again, 
the  soHd  line  shows  the  case  of  whitening  case  whereas  the  dotted  fine  shows  the  un-whitened 
case. 


Figure  3-17.  MUSIC  Algorithm 
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Note  that  the  results  are  not  as  good  as  in  the  previous  case.  However,  one  can  see 
that  after  several  iterations,  the  locations  of  the  targets  are  clearly  estimated  using  MUSIC  in 
conjunction  with  the  whitening  functional.  Again  as  in  the  above  case,  the  sohd  line  shows  the 
results  of  the  whitening  process  while  the  dotted  line  shows  the  results  of  the  un-whitened 
case. 

3.4.2.2  -  ESPRIT  Operator 

We  have  selected  the  same  scenario  as  previously  in  which  we  considered  an  array 
consisting  of  14  sensors  spaced  at  half  wavelength  and  3  targets  located  at  40°,  50°  and  70°, 
with  an  SNR  of  5  dB.  The  AR  process  was  assumed  to  be  of  order  2  with  parameters  1,  -0.9 
and  0.2.  The  following  are  the  results  of  50  runs. 


01=40” 

6,  =50” 

CD 

II 

o 

Mean 

41.1016 

51.3700 

70.6565 

Variance 

0.0018 

0.0026 

0.0046 

MSB 

60.6799 

93.8408 

21.5495 

Table  3-4.  ESPRIT  Estimates 


Note  again  that  the  bias  in  the  estimates  is  very  small.  However,  the  variances  of  these 
estimates  are  extremely  small.  Therefore,  one  can  see  the  advantages  of  this  technique  as 
compared  to  other  techniques  such  as  MUSIC.  Note  also  that  these  estimates  are  very  good 
because  of  the  fact  the  we  used  an  SVD  decomposition  on  the  two  sub-matrices  involved  in 
the  matrix  pencil.  Again  as  we  did  before,  we  can  visualize  these  results  by  looking  at  the 
angle  estimates  as  a  function  of  the  number  of  iterations  and  the  histograms  of  the  estimates. 


134 


Angle  =  50  Degrees  Angle  =  40  Degrees 
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Number  of  Occurences 


Figure  3-20.  Angle  =  70°  as  a  function  of  iteration  number 


136 


Number  of  Occurences  Number  of  Occurences 


Angle  =  50  Degrees 

Figure  3-22.  Angle  =  50°  (Histogram) 


Angle  =  70  Degrees 


Figure  3-23.  Angle  =  70°  (Histogram) 
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Even-though  the  angle  estimates  have  a  small  bias,  the  variances  of  these  estimates  are 
very  small. 

3.4.3  -  Signal  Plus  Clutter  Plus  Noise  With  Estimated  Whitenin2  Functional 

In  the  following  figures,  the  signal  was  assumed  to  be  of  the  form 

Y=AS  +  CN,+N2, 

where  Ni  and  N2  are  independent  white  noise  processes.  We  only  consider  the  case  where  we 
estimated  the  functional  and  apphed  it  to  the  AOA  algorithms  being  used.  To  achieve  this,  we 
first  estimated  the  AR  parameters  and  computed  the  AR  covariance  matrix  from  which  the 
whitening  functional  was  obtained.  The  following  figures  show  the  results  of  the  computer 
simulations.  Note  that  from  each  case,  we  first  present  the  values  of  the  AR  parameter 
estimates  used  in  the  estimation  of  the  angles  of  arrival,  then  the  corresponding  plots  or  tables. 
3.4.3.1  -  MUSIC  Operator 

In  this  section  we  present  the  results  of  the  MUSIC  algorithm. 
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AR  Parameters 


Figure  3-25.  MUSIC  Algorithm  with  p  =  0.6 


Number  of  Iterations 

Figure  3-26.  AR  Parameter  Estimates  (  p  =  0.7) 
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AR  Parameters  Spectrum  (dB) 


MUSIC  Algorithm  (Estimated  Whitening  Functional) 


Figure  3-27.  MUSIC  Algorithm  with  p  =  0.7 


Figure  3-28.  AR  Parameter  Estimates  (  p  =  0.6) 
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AR  Parameters  |  Spectrum  (dB) 


MUSIC  Algorithm  (Estimated  Whitening  Functional) 


Figure  3-29.  MUSIC  Algorithm  with  p  =  0.6 


Figure  3-30.  AR  Parameter  Estimates  ( o  =  0.7) 
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MUSIC  Algorithm  (Estimated  Whitening  Functional) 


Figure  3-31.  MUSIC  Algorithm  with  o  =  0.7 


Because  the  AR  parameters  have  been  estimated  using  the  technique  described  above 
as  well  as  the  corresponding  covariance  matrix  from  which  the  whitening  functional  was 
derived,  one  would  expect  the  results  to  be  slightly  different  than  those  of  the  known 
likelihood  functional.  To  our  surprise,  the  results  are  not  as  bad  as  one  would  expect  as  can 
be  seen  in  Figures  (3-24)  through  (3-31).  The  locations  of  the  signals  were  estimated  quite 
accurately,  except  for  small  biases.  This  clearly  shows  the  efficacy  of  the  proposed  technique. 
3.4.3.2  -  ESPRIT  Operator 

We  have  selected  the  same  scenario  as  previously  in  which  we  considered  an  array 
consisting  of  14  sensors  spaced  at  half  wavelength  and  3  targets  located  at  40°,  50°  and  70°, 
with  an  SNR  of  0  dB.  The  AR  process  was  assumed  to  be  of  order  2  with  parameters  1,  -0.9 
and  0.2.  The  following  are  the  results  of  50  runs. 
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Mean 


Variance 


MSE 


01=40“ 

02=50“ 

03=70“ 

39.8638 

49.9485 

69.9631 

0.2689 

0.1285 

0.2760 

0.9268 

0.1327 

0.0679 

Table  3-5.  ESPRIT  Estimates 


Note  again  that  the  bias  in  the  estimates  is  very  small.  However,  the  variances  of  these 
estimates  are  extremely  small.  Therefore,  one  can  see  the  advantages  of  this  technique  as 
compared  to  other  techniques.  Note  also  that  these  estimates  are  very  good  because  of  the 
fact  the  we  used  an  SVD  decomposition  on  the  two  sub-matrices  involved  in  the  matrix 
pencil.  Again  as  we  did  before,  we  can  visualize  these  results  by  looking  at  the  angle 
estimates  as  a  function  of  the  number  of  iterations  and  the  histograms  of  the  estimates. 
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Angle  =  50  Degrees 
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Number  of  Occurences 


LeCadre  Technique  (ESPRIT) 


39  40  41 

Angle  =  40  Degrees 


Figure  3-35.  Angle  =  40°  (Histogram 


Figure  3-37.  Angle  =  70°  (Histogram) 


Even-though  the  angle  estimates  have  a  small  bias,  the  variances  of  these  estimates  are 
very  small 

3.5  -  Effects  of  Correlation  in  the  AR  coefficients 

We  then  investigated  the  performance  of  the  technique  in  estimating  the  AR 
parameters.  It  was  observed  that  the  estimation  of  these  parameters  depends  greatly  on  the 
correlation  between  these  parameters.  For  example,  in  the  case  where  the  parameters  are 
uncorrelated  (  or  very  small  correlation),  the  algorithm  (minimization  of  the  whitened 
likelihood  functional  with  respect  to  the  AR  parameters)  performed  very  well  in  the  sense 
good  estimates  were  obtained  just  after  few  iteration  (10-15)  depending  on  the  scenario  used. 
However,  in  the  case  where  the  AR  parameters  were  strongly  (or  highly)  correlated,  the 
algorithm  did  not  perform  as  well,  since  over  100  iterations  are  needed  to  get  good  estimates 
of  these  parameters.  Furthermore,  the  estimates  are  biased. 
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Conclusions  and  Future  Recommendations 
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A  new  procedure  referred  to  as  A’SCAPE  has  been  introduced  in  the  first  part  of  this 
report.  The  goal  of  A’SCAPE  is  to  partition  a  scene  into  homogeneous  regions  and 
characterize  statistically  the  different  regions.  As  explained  in  Chapters  2  and  3,  two 
procedures  are  used  to  partition  the  scene.  First,  a  Mapping  procedure  is  used  to  separate 
between  regions  of  different  average  magnitude  levels.  This  procedure  is  a  new  adaptive 
thresholding  technique  which  proves  to  be  very  effective  at  separating  regions  even  when 
their  average  magnitude  levels  are  very  close  to  the  point  that  the  area  of  overlapping 
between  their  histograms  is  large.  Then,  when  all  regions  of  different  average  magnitude 
levels  have  been  separated,  a  Statistical  procedure  is  used  to  separate  between  sub-regions  of 
different  data  distributions  within  each  region  declared  as  homogeneous  by  the  Mapping 
procedure.  The  Statistical  procedure  is  based  on  the  Ozturk  algorithm  which  needs  only  100 
samples  to  properly  approximate  the  PDF  of  a  region.  In  the  Statistical  procedure,  Gaussian 
and  non-Gaussian  regions  are  formed  within  each  region  declared  as  homogeneous  by  the 
Mapping  procedure.  Outliers  are  detected  which  may  represent  small  patches  with  not 
enough  pixels  to  be  characterized  as  sub-regions.  Also,  contiguous  non-Gaussian  regions  that 
can  be  approximated  by  the  same  PDF  are  grouped  to  form  the  subpatches  within  each  region 
previously  declared  as  homogeneous  by  the  Mapping  procedure.  Furthermore,  PDFs  of  the 
non-Gaussian  regions  are  approximated. 

Though  using  two  procedures,  A’SCAPE  needs  four  stages  to  achieve  its  aim  at 
separating  the  different  contiguous  non-homogeneous  regions  that  may  exist  in  a  scene. 
These  stages  consist  of  a  Preprocessing  stage  where  classical  time-frequency  processing  of 
the  data  is  performed,  a  Mapping  procedure  stage,  a  Statistical  procedure  stage,  and  an 
Indexing  stage  which  assigns  a  set  of  descriptors  to  every  pixel  in  the  scene  under 
investigation.  When  A’SCAPE  is  followed  by  a  detection  stage,  all  information  is  available 
with  respect  to  which  detector  should  be  used  for  every  pixel  in  the  scene. 

Interactions  between  the  different  stages  of  A’SCAPE  are  controlled  through  an 
expert  system  shell,  IPUS.  This  is  important  because  the  goal  behind  A’SCAPE  is  to  monitor 
an  environment  which  is  unknown  a  priori.  Thus,  IPUS  is  needed  to  ensure  that  the  setting  of 
the  control  parameters  of  the  different  algorithms  in  A’SCAPE  is  well  suited  with  the  type  of 
data  being  analyzed.  If  not,  a  discrepancy  is  detected  and  IPUS  searches  for  the  distortion  that 
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may  have  caused  it.  This  results  in  the  choice  of  another  setting  of  the  control  parameters  and 
in  the  reprocessing  of  the  data. 

The  A’ SCAPE  procedure  is  illustrated  through  an  example  of  real  IR  data  collected 
over  a  lake  and  land  regions.  First,  the  preprocessing  stage  selects  those  pixels  in  the  scene 
that  result  in  uncorrelated  data  and  feeds  them  to  the  mapping  procedure.  Then,  the  mapping 
procedure  partitions  the  scene  into  three  main  patches;  lake,  land,  and  a  suhpatch  within  the 
land.  By  analyzing  the  statistics  of  the  three  regions  it  is  concluded  that  the  sub-patch  within 
the  land  may  be  a  body  of  water.  Following  that,  the  statistical  procedure  partitions  further 
the  lake  and  land  into  main  Gaussian  patches  and  a  total  of  fifteen  non-Gaussian  subpatches. 
Also,  outliers  representing  tiny  patches,  such  as  roads,  have  been  located  and  approximating 
PDFs  have  been  determined  for  the  non-Gaussian  subpatches. 

A  demo  package  built  in  Matlab  is  included  which  describes  in  detail  the  different 
stages  of  A’ SCAPE.  The  package  has  a  friendly  mouse  driven  graphical  user  interface  (GUI) 
and  consists  of  two  main  sections.  The  first  section  presents  the  detailed  steps  of  the  real  IR 
data  example.  The  second  section  is  subdivided  into  two  subsections  where  in  the  first  one  a 
set  of  examples  is  presented  which  illustrate  the  need  for  A’ SCAPE,  whereas  in  the  second  a 
description  is  given  for  every  stage  of  A’ SCAPE.  The  views  can  be  displayed  manually  or  in 
an  automated  way  as  a  slide-show.  A  movie  is  included  which  shows  the  steps  the  scene  goes 
through  during  the  mapping  procedure.  This  demo  package  is  among  the  deliverables. 

Work  under  investigation  and  future  work  include: 

(1)  Tailoring/tuning  of  the  edges  for  subpatches  detected  by  the  statistical  procedure 
which  is  limited  by  the  requirement  of  100  reference  pixels  for  the  Ozturk  algorithm. 
Note  that  when  the  subpatches  have  more  than  100  pixels  it  is  possible  to  do  more 
processing  to  fine  tune  the  edges. 

(2)  Study  of  the  performance  improvement  of  the  non-Gaussian  detectors  over  the 
Gaussian  detector  in  different  types  of  environments  and  under  different 
circumstances. 

(3)  development  of  more  expert  system  rules  to  enable  A'SCAPE  to  be  applied  to 
different  types  of  data  (e.g.,  radar,  IR,  sonar,  medical  imaging,  etc.). 
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(4)  Application  of  A'SCAPE  to  medical  imaging  (e.g.,  detection  of  tumors  in  lung),  to 
law  enforcement  (e.g.,  detection  of  concealed  weapons  imdemeath  clothing)  and  to 
other  areas. 

(5)  Design  of  a  detection  stage.  As  mentioned  before,  if  A’SCAPE  is  to  be  followed  by  a 
detection  stage,  the  type  of  detector  is  selected  according  to  the  approximating  PDF 
found  for  the  patch  under  investigation  and  parameters  for  the  sufficient  statistic  of 
the  detector  are  determined  from  the  pixel  descriptors  available  at  the  indexing  stage. 
Work  needs  to  be  done  to  test  the  detection  stage,  study  its  performance,  and  define 
expert  system  rules  to  overcome  any  possible  discrepancies. 

(6)  Incorporation  of  A’SCAPE  into  Rome  Laboratory  Space-Time  Adaptive  Processing 
(RL/STAP)  algorithms  to  enable  the  characterization  of  the  scene  under  investigation 
prior  to  any  detection.  As  mentioned  in  this  report,  when  A’SCAPE  is  done 
processing  the  scene  all  information  is  available  at  the  indexing  stage  for  every  test 
pixel  in  the  scene  with  respect  to  which  PDF  can  approximate  the  data  in  the  pixel, 
which  patch  does  it  belong  to,  which  pixels  can  be  chosen  as  reference  for  the  test 
pixel,  etc.  All  this  information  is  necessary  to  enable  proper  choice  of  the  detector 
type  and  its  corresponding  parameters  at  the  detection  stage. 

The  work  presented  in  Part  II  of  this  report  can  be  summarized  as  follows:  First,  we 
introduced  the  problem  of  detection  using  non-parametric  techniques  and  showed  the 
disadvantages  that  are  inherent  in  these  techniques.  Then  we  introduced  the  concept  of 
Parametric  detection  schemes  and  demonstrated  the  need  to  use  such  techniques  because  of 
the  number  of  degrees  of  freedom  they  add.  We  derived  the  beamforming  algorithm  based 
on  LeCadre’s  work  and  extended  to  include  both  clutter  and  additive  white  noise.  It  is  the 
parameterization  of  the  clutter  which  gives  us  the  driving  factor  to  include  the  added  degrees 
of  freedom  from  which  the  order  of  the  AR(MA)  model  used  is  estimated  and  then  estimate 
the  AR(MA)  coefficients  to  obtain  the  AR(MA)  covariance  matrix.  This  matrix  will  be  used 
to  “whiten”  the  covariance  matrix  of  the  received  vector.  Then  the  number  of  desired  signals 
present  in  the  environment  is  estimated  using  a  modified  version  of  the  so-called  MDL 
algorithm.  Once  these  parameters  have  been  estimated,  then  any  direction  finding  operator 
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could  be  used  to  determine  the  locations  of  the  signals.  Note  that  only  spatial  processing  has 
been  done  at  this  stage. 

In  the  simulations,  we  started  by  checking  the  convergenee  of  the  teehnique  and  its 
ability  to  estimate  correctly  the  AR(MA)  parameters.  We  found  that  the  teehnique  is  very 
efficient  in  obtaining  good  estimates  of  the  AR(MA)  parameters  with  few  iterations.  Note 
that  this  was  the  case  of  uncorrelated  AR(MA)  coefficients.  In  the  case  of  highly  correlated 
coeffieients,  as  many  as  500  iterations  are  needed  to  get  good  estimates  of  the  AR(MA) 
parameters.  Note  also  that  in  order  to  get  these  estimates,  the  order  of  the  AR(MA)  model 
has  to  be  known  a  priori.  In  most  cases,  it  is  estimated  from  the  data  set  along  with  the 
number  of  signals  present  in  the  environment.  We  actually  showed  that  an  extension  of  the 
MDL  algorithm  gave  the  exact  model  order  and  number  of  signals  present. 

We  then  tested  the  new  technique  in  location  estimation.  For  this,  we  used  two 
operators,  namely  the  MUSIC  operator  and  the  ESPRIT  operator.  The  reasons  for  ehoosing 
these  operators  is  purely  arbitrary.  Recall  that  the  first  operator  is  called  a  “search  procedure” 
operator  since  we  have  to  search  over  the  entire  array  manifold  to  obtain  the  DOA  estimates. 
This  may  be  very  computationally  extensive  and  the  amount  of  memory  needed  to  store  the 
results  may  be  very  high.  For  these  reasons,  we  chose  an  alternative  technique,  referred  to  as 
a  “non  search  procedure”  to  estimate  the  DOA’s  of  the  signals.  It  is  well  known  in  ESPRIT 
that  the  DOA  estimates  are  obtained  from  the  generalized  eigenvalues  of  the  a  matrix  pencil 
obtained  from  the  received  data. 

We  ran  both  algorithms  and  as  expected,  we  obtained  excellent  results  when  using  the 
“whitening  functional”.  Failure  to  do  so  results  in  poor  estimation  of  the  DOA’s.  We  have 
even  missed  some  targets  as  can  be  seen  from  the  results  presented  in  section  3. 

It  is  our  belief  that  the  newly  developed  technique  has  great  potential  in  detecting  the 
number  of  signals  as  well  as  determining  their  locations.  It  is  our  recommendation  that  it  be 
incorporated  into  the  Rome  Laboratory  Space-Time  Adaptive  Processing  (RLSTAP) 
capabilities. 
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In  the  future,  we  will  include  all  the  temporal  effects  to  the  technique  to  make  it  a  full 
space-time  adaptive  technique.  It  is  also  our  intent  to  include  the  effects  of  mutuals  coupling 
which  are  inherent  when  considering  an  array  of  sensors.  It  is  our  hope  and  belief  that 
compensation  technique  will  definitely  improve  the  quality  of  the  estimates  imder  study. 
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MISSION 

OF 

ROME  LABORATORY 


Mission.  The  mission  of  Rome  Laboratory  is  to  advance  the  science  and 
technologies  of  command,  control,  communications  and  intelligence  and  to 
transition  them  into  systems  to  meet  customer  needs.  To  achieve  this, 
Rome  Lab: 


a.  Conducts  vigorous  research,  development  and  test  programs  in  all 
applicable  technologies; 

b.  Transitions  technology  to  current  and  future  systems  to  improve 
operational  capability,  readiness,  and  supportabilHy; 

c.  Provides  a  full  range  of  technical  support  to  Air  Force  Materiel 
Command  product  centers  and  other  Air  Force  organizations; 

d.  Promotes  transfer  of  technology  to  the  private  sector; 

e.  Maintains  leading  edge  technological  expertise  in  the  areas  of 
surveillance,  communications,  command  and  control,  intelligence,  reliability 
science,  electro-magnetic  technology,  photonics,  signal  processing,  and 
computational  science. 

The  thrust  areas  of  technical  competence  include:  Surveillance, 
Communications,  Command  and  Control,  Intelligence,  Signal  Processing, 
Computer  Sdence  and  Technology,  Electromagnetic  Technology, 
Photonics  and  Reliability  Sciences. 


